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ABSTRACT 



The viscosity of the QGP is a presently hotly debated subject. Since its computation 
from first principles is difficult, it is desirable to try to extract it from experimental data. 
Viscous hydrodynamics provides a tool that can attack this problem and which may work 
in regions where ideal hydrodynamics begins to fail. 

This thesis focuses on viscous hydrodynamics for rclativistic heavy ion collisions. We 
first review the 2nd order viscous equations obtained from different approaches, and then 
report on the work of the Ohio State University group on setting up the equations for causal 
viscous hydrodynamics in 2+1 dimensions and solving them numerically for central and 
noncentral Cu+Cu and Au+Au collisions at RHIC energies and above. We discuss shear 
and bulk viscous effects on the hydrodynamic evolution of entropy density, temperature, 
collective flow, and flow anisotropies, and on the hadron multiplicity, single particle spectra 
and elliptic flow. Viscous entropy production and its influence on the centrality dependence 
of hadron multiplicities and the multiplicity scaling of eccentricity-scaled elliptic flow are 
studied in viscous hydrodynamics and compared with experimental data. The dynamical 
effects of using different versions of the Israel-Stewart second order formahsm for causal 
viscous fluid dynamics are discussed, resolving some of the apparent discrepancies between 
early results reported by different groups. Finally, we assess the present status of constraining 
the shear viscosity to entropy ratio of the hot and dense matter created at RHIC. 
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Chapter 1. Introduction 



1.1 The Quark Gluon Plasma and Relativistic Heavy Ion Collisions 

Shortly after the discovery of the asymptotic freedom of QCD, people realized that common nuclear 
matter, confining quarks and gluons within individual protons and neutrons, could transform into a new 
de-confinement phase (at high energy densities) - the quark gluon plasma (QGP) a form of matter 

that once existed in the very early universe or some variant of which possibly still exists in the inner core 
of a neutron star. It was theoretically conjectured that such extreme conditions can also be realized on 
earth through colliding two heavy nuclei with ultra- relativistic energies [3] , which transform a fraction of the 
kinetic energies of the two nuclei into heating the QCD vacuum within an extremely small volume. However, 
it actually took more than 25 years of efforts to reach the threshold for the QGP phase transition, from the 
first relativistic heavy ion program at the Bevalac at LBL (with beam energies of up to IGeV/nucleon) in 
the early 1980s to the AGS at BNL (with center of mass energies per nucleon pair y/s ^ 5GeV) and the 
SPS at CERN (with - 17GeV) in the late 1980s 0, to the RHIC at BNL (with = 200GeV), starting 
in 2000 [51[71[51[g|. and to the inpending heavy ion program at the LHC at CERN (with = 5.5TeV), which 
will probably start towards the end of next year [TO] . 

At AGS and SPS energies and below, there was no unambiguous evidence for QGP formation, although 
a number of signals found at the SPS strongly suggested the formation of a "new state of matter" pl^H 
Only after the beginning of the RHIC program during the summer of 2000, more and more evidence showed 
that the QGP had been discovered [6l[7l[8l[9l fT2l[T3l[T4] . Before discussing these QGP signatures, we will first 
explain the concept of the QGP and how its features have become richer during the past few years. 



1.1.1 Prom weakly coupled QGP to strongly coupled QGP 

The quark gluon plasma is defined as a thermalized state of quarks and gluons without color confinement. 
It was originally conceived as a weakly coupled gas, motivated by the asymptotic freedom of QCD at high 
energies. Using the statistics of relativistic massless fermions and bosons, one obtains the equation of state 
(EOS) for a free massless QGP gas [15] : 
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Here, gg + ^{gq + gq) is the total degeneracy of the QGP. For gluons with 8 colors and 2 polarizations, 
gg — % X 2. For quarks with 3 colors, 2 spins and Nf flavors, = 3 x 2 x Nf. The 7/8 factor comes from 
Fermi-Direc statistics. The above equations give the ideal EOS for a non-interaction massless QGP, p = e/3. 
If one assumes that the phase transition happens at e « IGeV/fm'^, roughly seven times that of normal 
nuclear matter, then one finds a critical temperature of w 160 MeV for a massless 2 flavor (u, d) QGP. 

In principle, the free energy and the EOS of the weakly coupled QGP can be perturbatively calculated 
order by order in thermal QCD. However, the QCD free energy shows bad convergence of the perturbative 
expansion in powers of the strong coupling constant g |16pi7j . To solve this problem, one needs to reorganize 
the perturbation theory. This led to the recent developments of hard thermal loop (HTL) perturbation 
theory, the so called 0-derivable approach, dimensionally reduced screened perturbation theory, and others 



^ The SPS measurements indicate that the plasma was already created there, but was very short lived. 
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Figure 1.1: Energy density and three times the pressure as a function of temperature, calculated from lattice 
QCD with asqtad and p4 action pT] , 

(see the review article fTB| for details) , all of which yield much improved convergence for the QCD free energy 
and an improved EOS at high temperature. Still, the calculation needs non-perturbative input, due to the 
existence of a non-perturbative magnetic mass scale nimag ~ 9^T, which enters the EOS at order [roil20j . 
On the other hand, even the most sophisticated weakly coupled QCD methods fail near the phase transition 
where non-perturbative effects become dominant. 

Lattice QCD offers a non-perturbative approach to study the QCD properties at finite temperature. Re- 
cent developments in this field have made available quite precise lattice calculations of the EOS with almost 
physical quark masses [HJI^^]. Fig 11.11 shows the Lattice EOS (at zero chemical potential) from Bazavov 
et al., using two different improved staggered fcrmion actions (the asqtad and p4 action, respectively) [21j . 
They found that both deconfinement and chiral symmetry restoration happen in a narrow temperature re- 
gion: 185 MeV <Tc< 195 Mev|, which is indicated by the narrow band in Fig ll.ll At higher temperature 
T > 300 MeV, one finds that the lattice EOS approaches 85-90% of the Stefan-Boltzmann limit for an ideal 
non-interacting QGP gas. But this does not necessarily mean that the QGP is already weakly coupled there. 
From the EOS alone, lattice QCD simulations can not distinguish a weakly coupled QGP from a strongly 
coupled QGP. Indeed, the concept of a strongly coupled QGP, which behaves like an almost perfect liquid 
with very low viscosity |24[12[13p25j , came from the strong collective flow observed in RHIC experiments and 
its very successful descriptions by ideal hydrodynamics [SSIEZ]- However, why the QGP is strongly coupled 
at RHIC energies is still a theoretical challengq^l, and it is also unknown whether the QGP created in the fu- 
ture LHC experiments will continue to be a strongly coupled liquid or behave more like a weakly coupled gas. 

1.1.2 Stages of a heavy ion collision and theoretical tools 
Stages of a heavy ion collision 

In this subsection, we introduce different stages of a heavy ion collision and the main theoretical tools to 
describe or simulate these stages. After a discussion of these theoretical tools, we will list several experimental 
probes and signatures to detect the transient QGP phase. 

^ Other groups found a lower transition temperature 155 MeV < Tc < 174 MeV 

^ Some models trying to explain the strongly coupled nature of the QGP at RHIC energies can be found in Refs. \28l\29f . 
Other authors explore the possibility that the strongly coupled dynamics can be obtained through weak coupling expansions using 
HTL resummation fj3(M31}/ . 
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Figure 1.2: Different stages for relativistic heavy ion collisions [33l 



The different stages for an ultra-relativistic heavy ion collision {y/s > 100 GeV) are schematically illus- 
trated in Fig 11.21 and 11.31 The initial stage before the collision is followed after impact by a pre-equilibrium 
stage, an expanding QGP and hadron resonance gas (HRG) stage and a final freeze-out/decoupling stage. 
In more detail, these stages are characterized as follows (32; : 

1. Initial stage at r < 0: Two Lorentz contracted heavy nuclei approach each other with more than 99.9% 
of the speed of light. At sufficiently high collision energy (possibly reached at RHIC), this initial stage 
can be described by dense gluon walls known as the Color Glass Condensate (CGC) (this is further 
explained in the theoretical tools below). 

2. Pre-equilibrium stage and thermalization: The energetic collision of the two heavy nuclei excites the 
QCD vacuum and produces a dense pre-equilibrium matter consisting of quarks, anti-quarks and 
gluons. It takes around 1 fm/c for the pre-equilibrium bulk matter to achieve local thermalization and 
form the quark-gluon plasma. In this very early pre-equilibrium stage, the primary collisions between 
fast partons inside the colliding nuclei also generate "hard probes" with either large mass or large 
transverse momentum, such as heavy quark pairs (cc and bb) , pre-equilibrium real or virtual photons, 
and very energetic quarks and gluons with large transverse momentum (from which jets are formed 
after hadronization) . 

3. QGP expansion and hadronization: After thermalization, the quark-gluon plasma (QGP), driven by 
thermal pressure gradients, expands and cools down very quickly. After reaching the critical temper- 
ature Tc ~ 170 MeV, it hadronizes and turns into hadronic matter, consisting of a mixture of stable 
and unstable hadrons and hadron resonances. Hadronization happens continuously at the edge of the 




Figure 1.3: The space-time picture of a heavy ion collision [51] . 
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QGP fireball during the whole QGP expansion period. In central Au+Au collisions at RHIC, it takes 
around 10 fm/c for the QGP fireball to expand and fully convert to hadronic matter. 

4. Hadronic expansion and decoupling: The hadronic matter continues to expand until the system be- 
comes very dilute. Then individual hadrons decouple from the system (kinetic freeze-out) and free- 
stream to the detector. Like the QGP hadronization process, the hadronic decoupling happens con- 
tinuously at the edge of the fireball, where the density is low. After complete hadronization, it takes 
another 5—10 fm/c for the hadronic matter to completely freeze-out. 

Theoretical tools for different stages of a heavy collisions 

There does not exists a unique theoretical tool to describe the whole heavy ion collision process from the 
very beginning till the end. The different energy and time scales during different collision stages imply dra- 
matic changes of the effective physical degrees of freedom and their interactions, which thus require different 
tools for their description. A complete description requires matching these tools to each other, generating 
so-called hybrid approaches. For example, the results from the Color Glass Condensate (CGC) theory can 
be used as initial conditions for dynamical evolution models such as hydrodynamics or the parton cascade. 
Hydrodynamics or parton cascade models are then connected with a hadron cascade model for a description 
of the late hadronic stage. The hadron cascade model is equipped with an afterburner to generate quantum 
statistical or final-state interaction induced 2-particle correlations. Fig 11.41 shows the different theoretical 
approaches and their ranges of applicability. Although plotted by S. Bass as early as 2001, it still very nicely 
illustrates our present understanding, in spite of several new theoretical developments in the past yearfl 
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Figure 1.4: Theory tools for RHIC and their range of applicability [39]. Solid bands denotes the safe 
range to apply respective theoretical tools, while the dashed and dotted bands refer to the region where the 
approach is still applied, but may be questionable or unsafe. GYM: Classical Yang-Mills Theory [40ll41j. 
LGT: Lattice Gauge Transport [42, 43j, PCM: Parton Cascade Model gj, NFD: Nuclear Fluid Dynamics 
or Hydrodynamics [45[[26ll27] . String and Hadronic TM: String and Hadronic Transport Model [46ll47l[48] . 



* These new developments include new ideas to understand thermalization f35ll 30f as well as first numerical implementations 
of hybrid approaches, such as hydrodynamics +hadron cascade 
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Color Glass Condensate: At very high colhsion energies, particle production at mid rapidity probes 
the nuclear structure functions in the small- a; regime (a; is the fraction between parton momentum and beam 
momentum). It is well known that the gluon distribution function xG{x,Q^) increases dramatically with 
decreasing of x, for large enough probe resolution Q^. When the gluon density becomes high enough, two 
gluons start to recombine to one. This leads to gluon saturation below some momentum scale Q^. At this 
momentum scale, each gluon mode has macroscopic occupation number ~ which is why this state has 
been called a condensate. These small- a; gluons are generated by radiation corrections from gluons at large 
X, whose natural evolution time scale is Lorentz dilated. This time dilation is transferred to the small- a; 
degrees of freedoms, making them evolve very slowly compared to other natural time scales: their behavior 
is "glassy" . Considering these factors, together with the color nature of gluons, this initial stage has been 
named Color Glass Condensate (CGC) [40l|4TJ|49]. The production of the pre-equilibrium secondary soft 
gluons is related to breaking the phase coherence in this initial CGC wave function [50] . 

Parton cascade model: After the initial parton production, one needs to describe the pre-equilibrium 
stage and its thermalization. For dilute systems with incoherent parton configurations, the classical motion 
of on-shell partons can be described by the Parton Cascade Model (PCM) [44l[5Tl[52], which solves the 
Boltzmann equation with a leading order PQCD collision term. Initially, the PCM was assumed to work for 
the pre-equilibrium stage, the subsequent thermalization period and the succeeding QCP expansion stage. 
However, the typical thermalization time obtained from PCM with 2-body {gg gg) PQCD scattering cross 
sections is of the order of 5 fm/c [52], which is too long for the fast thermalization (< 1 fm/c) required by 
RHIC data [2ll[53]. The recent development of a PCM with 2^3 {gg ^ ggg) processes leads to faster 
thermalization and indeed reproduces the large observed elliptic flow f54','55]. However, it is applied to a 
dense system, which creates tension with the assumptions under which the Boltzmann equation is valid. In 
a very dense system (such as the very early pre-equilibrium stage at RHIC and LHC energies) the partons 
scatter so frequently that they are no longer on shell. To treat partons with off-shell energies requires the use 
of quantum transport theory. The theoretical framework for quark-gluon quantum transport was developed 
by Heinz more than 20 years ago [5Sl[571[SHl[5n] , but its numerical demands are exorbitant, and it has not 
yet been implemented numerically. 

Thermalization: The formation of the quark gluon plasma requires two aspects: local equilibrium 
(thermalization) and local momentum isotropy. Generally, it is believed that the system achieves momen- 
tum isotropy before completing thermalization, which is assumed to happen at a time scale of ~ 0.6 fm/c 
at RHIC energies as indicated by the validity of ideal hydrodynamic simulations [24l[53]. Elucidating the 
mechanism for fast isotropization and thermalization has been a theoretical challenge. Although much 
progress has been made during the past years, the thermalization mechanism is still not fully understood. 
Recently, people realized that the Weibel instability, which is known from traditional electromagnetic plas- 
mas |6^, might help to understand the fast isotropization and thermalization of the QGP [6n[62lf35l l36j . 
The presence of a Chromo- Weibel instability in the early parton plasma, which is characterized by strong 
momentum anisotropics, has been confirmed by numerical simulations in 1 -Hi-dimensional |35 1 l36 ( f63] and 
3-|-l-dimensional [64,, 65 hard thermal loop effective theory and in a 3-|-l-dimcnsional expanding Glasmad 
[66l[67] . However, non-Abelian saturation effects temper the exponential growth of the Chromo- Weibel in- 
stability, leading to significantly larger time scales even for isotropization than those required for the fast 
thermalization approximately observed at RHIC [64ll66j. 

Hydrodynamics: If thermalization is achieved and can be locally maintained during the subsequent 
expansion, the further evolution of the QGP and hadronic matter can be described by hydrodynamics [45] . 
Hydrodynamics is a macroscopic approach which describes the system by macroscopic variables, such as local 
energy density, pressure, temperature and flow velocity. It requires knowledge of the equation of state, which 

^ Glasma is the name introduced for initially highly coherent gluon matter that makes the transition from the Color Glass 
Condensate to the Quark Gluon Plasma iSOf . 
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gives a relation between pressure, energy and baryon density, but no detailed knowledge of the microscopic 
dynamics. The simplest version is ideal hydrodynamics |26[|27j . which totally neglects viscous effects and 
assumes that local equilibrium is always perfectly maintained during the fireball expansion. Microscopically, 
this requires that the microscopic scattering time is very much shorter than the macroscopic expansion (evo- 
lution) time and that the mean free path is much smaller than the system size. If this is not satisfied, viscous 
effects come in, and one can apply viscous hydrodynamics as long as the deviation from local equilibrium 
remains small |68[I69| . If the system is far away from equilibrium, one has to switch to a kinetic theory 
approach, such as parton [44j or hadron [47] cascade models. 

Hadron cascade model and hybrid approaches: The hadron cascade model [46ll47ll48j . which solves 
the Boltzmann equation for a variety of hadron species with flavor-dependent cross-sections, is a successful 
tool to describe the hadronic matter created at AGS and SPS energies. At these collision energies, the hadron 
cascade model is initialized by a superposition of hadrons and hadronic strings, produced in the primary 
nucleon-nucleon collisions. At RHIC energies and above, hybrid approaches that combine a parton cascade 
model or hydrodynamics with a hadron cascade, provide a "unified" description of the evolution of the QGP 
and the succeeding hadronic matter. However, some caution must be taken for the transition between the 
models. Parton -I- hadron cascade hybrids must deal with the problem of converting partons to hadrons 
without violating the second law of thermodynamics (i.e. without losing entropy), and they have difficulties 
to incorporate the change in the structure of the QCD vacuum during the phase transitioro. Hydrodynamics 
-I- hadron cascade hybrids [37l[38] can more easily accommodate these, by employing a realistic EOS from 
lattice QCD. One generally stops hydrodynamics at a switching temperature slightly below Tc, converting 
the fluid to hadrons using a Cooper-Frye prescription (see Chap. 3.4) for phase-space distributions to gener- 
ate (via Monte-Carlo) initial momentum and spectra profiles for the hadron cascade simulations. However, 
this procedure can not deal with a potential feed-back from cascade hadrons to the hydrodynamics fluid 
along space-like parts of the matching freeze-out hypersurface. 

1.1.3 QGP signatures 

Possible signals and probes for the quark-gluon plasma have been investigated for around 30 years since 
the birth of the field. Such signatures include: collective flow [7^73ll26ll?7] . strangeness enhancement [TiUTS] . 
charmonium suppression [751177] . thermal photon and dilepton emission [751175] . jet quenching [5Dll5n[5^[55] . 
critical fluctuations [84], and others. Some of the predicted signature were already found in earlier heavy ion 
experiments at AGS and SPS energies [11]. However, none of these signatures allow individually to prove 
QGP formation, as they are contaminated by the dynamical evolution of the fireball through various stages, 
usually from the very early pre-equilibrium stage through (perhaps) a QGP phase to the late hadronic stage. 
The combination of three observations at RHIC, that finally convinced the community that the QGP has 
been successfully created, were the measurements of strong anisotropic collective flow, valence quark number 
scaling of the elliptic flow V2, and jet quenching [85l[86]. This led to the announcement in 2005 that the 
QGP had been discovered at RHIC [12l[13l|85l|87] . 

Collective flow : 

The hadron momentum spectra, their angular distribution (flow patterns) and the particle yield ratios are 
primary observables for the bulk medium created in heavy ion collisions. One of the main discoveries of RHIC 
is that the medium displays strong collective dynamics [88118911901191] . which, for the first time in the history of 
particle and nuclear physics, could be quantitatively well described by ideal hydrodynamics[l|. In the language 

^ Some development try to solve this problem can be found in f70[|7-?y 

^Although large collective flow was also observed in heavy ion collisions at AGS and SPS energies, the corresponding data 
can not be quantitatively described by hydrodynamics models \92h 
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Figure 1.5: Elliptic flow V2 for different hadron species [94l[90l[95], plotted as a function of transverse 
momentum, compared with ideal hydrodynamics predictions |961l97j. 



of hydrodynamics, the collective flow is driven by pressure gradients, thus providing access to the equation 
of state (EOS) of the medium. Whereas the azimuthally averaged radial flow receives contributions from all 
expansion stages, the anisotropic "elliptic" flow seen in non-central collisions is generated mostly during the 
hot early stage, and thus provides the information about the QGP phase, namely its thermalization and its 

EOS [Miniinsi- 

The elliptic flow V2 is defined as the 2nd Fourier coefficient of the azimuthal distribution of hadron 
spectra: 

r dchr, COS(20„) 

V2{pTj = (cos(2</>p)) = , (1.2) 

J rp dypTdpTdfpp 

where ^.ypTdpT d4> ^® angular distribution of the transverse momentum (jpt) dependent spectra, and 
y = \ h\[{E +pl)/{E — pl)] is the rapidity of the particles. 

Fig 11.51 compares the experimental elliptic flow data V2{pt) with ideal hydrodynamic predictions |94[|90[ 
[951196] . For px < 1-5 GeV, where most (more than 98%) of the particles are produced, ideal hydrodynamics 
shows excellent agreement with the experimental data and correctly predicts the observed splitting for 
different hadron species. This strongly indicates that the bulk of the matter is strongly coupled and behaves 
like an almost perfect fluid [85] . 

For a successful description of the RHIC data, especially for the elliptic flow, ideal hydrodynamic requires 
a fast thermalization of the system, which must happen on a time scale of about 0.6 fm/c |24[l53j. Around 
that time, the early matter has a peak temperature of ~ 350 MeV - about twice the QGP phase transition 
temperature. Although this gives indirect evidence for QGP formation, the success of ideal hydrodynamics 
is primarily evidence for the formation of a thermalized new form of matter at T ~ 2Tc |24 1 I12| , [T3] . 

Quark degrees of freedom and partonic collectivity: 

The success of the statistical model in analyzing of the measured hadron abundances [§51112] gives ad- 
ditional evidence for the QGP phase transition. Particle ratios including a large number of hadron species 
are well described by a thermal model with just two parameters, T and /i, yielding a chemical freeze-out 
temperature Tch — 160 — 170 MeV, which is approximately equal to the QCD phase transition temperature 
as determined by lattice QCD simulations. This strongly indicates that hadrons are born during the QGP to 
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Figure 1.6: left: The transverse kinetic energy KEt dependence of V2 for various hadron species. Mesons 
and baryons fail to different universal curves respectively, right: V2 and KEt are scaled by the number of 
valence quarks. All of the hadron species fail to the same universal curve: the differential elliptical flow per 
quark. [TOU] 



hadron phase transition, and that flavor changing interactions (from inelastic collisions) quickly cease right 
after the phase transition [32] . 

Direct evidence for quark degrees of freedom in the newly formed matter can be extracted from detailed 
measurements of the p^-differential elliptic flow V2{pt) for a large variety of mesons and baryons. This 
observable, shown in Fig. 11.51 shows a characteristic splitting between mesons and baryons at intermediate 
transverse momentum px = 2—5 GeV. Even though this is above the pr range where hydrodynamics is valid, 
the underlying collective flow of the fireball still affects hadron production at this px [32|. After replacing 
the transverse momentum pT by the transverse kinetic energy KEt = ^/p^^+m^ — m, the mass splitting at 
Pt < 1.5 GeV disappears. Fig. 11.61 left shows, as a result of this procedure, two universal V2 scaling curves 
for baryons and mesons, respectively [lOOj . Baryons, which contain three valance quarks, show stronger V2 at 
intermediate transverse momentum than mesons, containing only two valance quarks. This phenomenon can 
be well explained by the quark recombination model [101.102 jll03pi04j . in which collectively flowing baryons 
and mesons are generated by the coalescence of quarks that collectively flow with the medium. According to 
the recombination model, the baryon and meson elliptic flow coefficients are expressed in terms of the quark 
elliptic flow in the following way 1105^ : 

vi''\pT) = 2vi'\pT/2)- vi''\pT) = 3vi''\pT/3). (1.3) 

As the right panel in Fig. ll.6l shows. a corresponding rescaling of both V2 and KEt by the number of valence 
quark (2 for mesons and 3 for baryons) leads to a universal V2 scaling curve for all hadron species. In short, 
the universal valence quark number scaling V2 suggests that the collectively flowing matter directly involves 
quarks, and that the quark collective properties are transferred to those of hadrons by quark recombination. 

Jet quenching: 

At the very beginning of a heavy ion collision, hard scatterings of incoming quarks and gluons every 
now and then create a pair of energetic fast partons with large transverse momentum. Each of the two fast 
partons will finally fragment into a spray of hadrons, forming what is called a jet. The rates for such hard 
process are small but grow rapidly with increasing collision energy. At RHIC energies, the fast partons for 
the first time became sufhciently abundant as useful probes for the hot medium. Their production rates are 
well understood, both experimentally from pp collision at the same energy or theoretically via PQCD. What 



8 



PHENIX Au+Au (central collisions): 




2 4 6 8 10 12 14 II 

Pj (GeV/c) 



Figure 1.7: The nuclear modification factor factor Raa as a function of transverse momentum px for direct 
photons 7, as well as 7r° and rj mesons in central Au + Au Collisions |106| . Raa is defined as the ratio of 
the cross section per nucleon-nucleon collision measured in a heavy ion collision divided by the cross section 
measured in p + p collisions. Raa = 1 if the heavy collisions can be viewed as a simple superposition of 
p + p collisions. 



makes them useful for heavy ion collisions is that, before hadronization to jets, they have to travel through 
the hot fireball medium formed in the collision, which modifies their initial energy and momentum. 

One of the most exciting results from RHIC, right after the discovery of strong elliptic flow, was the 
observation of a strong suppression of high-pr hadrons in central Au+Au collision, compared with the 
scaled results from the pp collisions [1091111011111] . The experimentally observed suppression by a factor 
4-5 agrees qualitatively with theoretical predictions for jet quenching [5Dllll2l[5^[55j . which argued that the 
QGP formation could lead to large energy loss of fast partons by collision induced gluon radiation and thus 
suppresses the production of high pT hadrons fragmented from such partons. The discovery of jet quenching 
is illustrated in Fig 11.71 The strong suppression of pions and rj mesons is distinctly in contrast to that of 
direct photons, which show no suppression since they interact only electrodynamically and thus directly 
escape from the medium without further interactions. 

Angular correlations between a high-pr leading (trigger) hadron with other energetic hadrons provide 
additional strong support for the picture of significant parton energy loss in the QGP medium 81j. Energy 
momentum conservation requires that fast partons are always generated in pairs, moving along opposite 
directions in the pair center of mass frame. This leads to back-to-back correlations between the resulting 
jets, as shown by two peaks in the azimuthal angle correlation separated by 180°. Such a back-to-back cor- 
relation is clearly seen in p-l-p and d-|-Au coUisions (see left panel of Fig. 11.81 |107[|113j ). In central Au-|-Au 
collisions, however, one sees only one such peak in the direction of fast trigger hadron (blue stars in the 
left panel of Fig. II. 8p . This can be understood if one assumes that the energetic parton pair is created 
near the surface of fireball: the near-side outgoing parton quickly escapes from the medium and fragments 
into the leading hadron and other softer hadrons correlated in angle with the leading hadron, while its 
inward-traveling partner at 180*^ loses most of its energy through interactions with the QGP medium and 
no longer contributes to energetic hadrons in the away-side (recoiling) direction [81j . The energy carried by 
the away-side fast parton is deposited in the medium, which leads to an enhancement of soft hadron pro- 
duction in the away-side hemisphere, as shown in the right panel of Fig ll.8l [108j . Soon, people realized that 
the much broadened away-side correlations may be related to a possible collective "hydrodynamic" response 
of the medium to the energy and momentum deposition from the fast parton, called Mach Cone |114[115|lll6j . 
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Figure 1.8: Azimuthal angular correlations between high momentum hadrons in p + p, d + Au, and central 
Au + Au collisions. In both cases, the trigge particle of the pair has high-momentum {pr > 4GeV). Left: 
Associated particles recoiling with high momenta {pr > 2GeV) exhibit strong suppression in Au + Au | 107j . 
Right: Associated particles with low recoil momenta {pr > O.lSGeV) are strongly enhanced in Au + Au 
collisions [108] . 

1.1.4 "RHIC scientists serve up the perfect liquid " 

The 2007 NSAC Long Range Plane for Nuclear Physics States that "The experiments performed at 
RHIC since it began operation in 2000 have provided spectacular evidence that the QGP does exist" [55] . 
but its properties are quite different from the earlier expectation based on PQGD, which had suggested that 
the QGP should behave like a dilute gas. In fact, the strong collective flow together with the very good 
description from ideal hydrodynamics indicates that the quark gluon plasma created at RHIC is strongly 
coupled and behaves like a nearly perfect liquid with very low viscosity [8_5j. 

The discovery of the strongly coupled QGP is intellectually exciting since it surprisingly connects super- 
string theory with relativistic heavy ion experiments. Mapping strongly coupled quantum field theories to 
weakly coupled gravity by the AdS/CFT correspondence |117[|118| . string theorists have developed tools 
for gaining (at least) qualitative insights into the strongly coupled QCD-like systems |119[|120| where tra- 
ditional PQCD methods fail to apply. Meanwhile, insights from other strongly coupled systems may help 
us to understand the nature of the strongly coupled QGP. An example comes from strongly interacting 
fermion systems, created in optical traps at extremely low temperatures, which show similar hydrodynamic 
behavior |12HI122] . The advantage of such cold atom experiments lies in the tunable coupling strength, via 
Feshbach resonances, which allows to switch between strongly and weakly coupled fermion systems and thus 
may help us to understand the transition between strongly and weakly coupled QGP. 

With these exciting developments in the past years, the next phase of the RHIC physics program and 
of the incoming heavy ion program at the LHC will focus on detailed investigations of the QGP, "both to 
quantify its properties and to understand precisely how they emerge from the fundamental properties of 
QCD" [85^. Fundamental questions that need to be addressed include (see Ref [123] for details): 

• What is the mechanism of the unexpectedly fast thermal equilibrium? 

• What is the initial temperature and thermal evolution of the produced matter? 

• What is the energy density and equation of state of the medium ? 

• What is the viscosity of the produced matter? 

• Is there direct evidence for deconfinement, color screening, and a partonic nature of the hot dense 
medium? What is the screening length? 
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• Is the chiral symmetry restored by QCD? 

• How does the new form of matter hadronize at the phase transition ? 

In this thesis, I wiU concentrate on one of the above questions, the viscosity of the QGP. To answer it 
requires continuous progresses on both theoretical and experimental sides, especially the development of vis- 
cous hydrodynamics and future high statistics flow measurements for a variety of identified hadrons species. 
This thesis will focus on viscous hydrodynamics for relativistic heavy collisions^. The transport coefficients 
(such as shear viscosity and bulk viscosity) are free parameters in viscous hydrodynamic calculations. The 
hope is that by tuning these parameters for best fits of a sufficiently large set of sensitive experimental 
observables, the QGP viscosity can be extracted phenomenologically from experimental data. This requires 
not only the development of a practical and accurate viscous hydrodynamics code, but a careful investiga- 
tion of other ingredients (EOS, initial and final conditions, etc.), which may affect the viscosity-sensitive 
observables. Some of these issues will be addressed in this thesis. (Of course, the QGP transport coefficients 
are not just simple numbers, but depend on the thermodynamics properties of the fireball, which evolve with 
time). We should therefore use, as much theoretical knowledge as available, to constrain the the temperature 
dependence of these variables, at least at the qualitative level. In the rest of this chapter we will briefly 
review the present status of theoretical understanding and knowledge of the transport coefficients for QCD 
matter and other strongly coupled systems. 



1.2 Transport coefRcients for non-relativistic fluids 

There are several transport coefficients that characterize the internal "friction" in a fluid. For example, 
the shear viscosity (dynamic viscosity) rj measures the fluid's resistance to flow, the bulk viscosity (volume 
viscosity) C measures the fluid's resistance to expansion, and the thermal conductivity (heat conductivity) A 
measures the fluid's ability to conduct heat. In this section, we will discuss shear and bulk viscosity in some 
detail, since they will be used in later parts of this thesis. 

Shear viscosity 

Classically, the shear viscosity is defined as the ratio between the friction force F per area and the 
transverse flow gradient VyVx, 

F 

j^rjVyVx. (1.4) 

Microscopically, shear viscosity is associated with momentum transfer between particles in different fluid 
cells. For a dilute gas of non-relativistic particles, the shear viscosity rj and shear viscosity to entropy density 
ratio rj/s are estimated as [147] : 

_ mv T] mv 

~ nmw/infp ~ , -~ , (1-5) 

a s an 

where n is the particle density, m is the particle mass, v is the mean velocity, /,„fp = (where a is the 
transport cross section) is the mean free path, and the entropy density s is proportional to the particle 
density s ^ n. The shear viscosity of a non-relativistic dilute gas increases with temperature 77 ~ V mT 
(taking v ^ y^T/m) and is approximately independent of densitjj^. The 1/n density dependence of 77/s 

® This thesis mainly covers work done and results obtained by the auth or \ 1 24f 1 251 1 26}\ 1 2Tf 1 28^12^ 1 30\l . Results from other 
groups published over the same time period can be found in Ref f jgjU J38U jg!^ jj^pgliTggy JgypflfTigpT^^ 
We uiill refer to these papers for comparison uihen needed or useful. 

® This estimate agrees qualitatively with results derived from the Chapman- Enskog expansion \14'Ti - 
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translates into a pressure dependence of 77/s, as shown in the left panel of Fig ll.9l Using the ideal gas EOS 
p = nkgT, one finds that, for fixed pressure, rj/s increases with temperature approximately as 

r)/s^T^/^ (1.6) 

in the high temperature gas phase. 

For a liquid, the momentum transfer between different fluid cells involves the motion of voids. The 
density of voids decreases with temperature, which leads to an increase of shear viscosity due to the growth 
of mean free path. Detailed calculation [148j shows that 

~ /ine^/C^-^) ^ ^ /,e^/(fc-^), (I.7) 
s 

where E is an activation energy and h is the Planck constant, and again s ~ n. This result shows that for a 
liquid rj/s increases with decreasing temperature. 

Taking the liquid phase and gas phase together, Eq. (|1.6|) and Eq. (|1.7p imply that rj/s is likely to reach a 
minimum during the liquid-gas phase transition. This is confirmed by experimental results shown in Fig 11.91 
which plots in the right panel r]/s as a function of temperature for three different fluids (helium, nitrogen 
and water) [146j . and for water at different pressures in the left panel [145j . In all of these cases, rj/s reaches 
a minimum near the phase transition. 

Bulk viscosity 

In non-relativistic fluid dynamics, the bulk viscosity C is defined as a combination of shear viscosity 
77 and volume viscosit}0 v. = v + ^rj. Bulk viscosity vanishes for a scale invariant system. For a 
dilute monatomic gas, experimental data and kinetic theory show that bulk viscosity is vanishingly small 
or zero [149] . The bulk viscosity of a dilute diatomic gas, however, has an appreciable value due to the 
exchange of energy between translational and rotational degrees of freedom during collisions [150] . In a 
dense gas, bulk viscosity is associated with the internal friction force arising from the change of volume at 
constant shape [151| . In a liquid, bulk viscosity is related to the rearrangement of molecules during acoustic 
compression and rarefaction [151j . 

Experimentally, bulk viscosity is generally determined from measuring the sound absorption coefficient. 
In contrast to shear viscosity, there are no comprehensive bulk viscosity data sets for many varieties of fluids 
over a wide range of temperature and density, and the existing bulk viscosity data have large error bars. 

The volume viscosity (also known as "second viscosity" ) v measures the resistance to expansion; it determines the dynamics 
of a compressible fluid. 
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The phase transition behavior of the bulk viscosity of spherical molecules was studied by Meier, Laesecke 
and Kabelac using molecular-dynamics simulations with Lennard-Jones potentials. They found that the 
bulk viscosity shows a peak in the vicinity of the gas-liquid phase transition |152j (see also Fig. 6 and related 
comments in Ref. |153j ). 

In short, the shear (bulk) viscosity to entropy ratio reaches a minimum (maximum) near phase transitions 
for common non-relativistic fluids. As we will see in Sec. II. 4i the same holds for relativistic QCD matter. 
This is not a rigorous statement from first principles calculations, but appears to be supported by specific 
examples of both experimental data and theoretical results. 



1.3 Transport coefRcients for relativistic fluids: weak coupling vs. strong cou- 
pling 

Relativistic hydrodynamics can be viewed as an effective theory for quantum field systems at large 
distance and time scales. The dynamics of long wavelength, low frequency fluctuations are characterized by 
transport coefficients (shear viscosity 77, bulk viscosity C and heat conductivity A, electric conductivity a etc.). 
For a relativistic fluid, the transport coefficients define the leading order corrections of the energy momentum 
tensor T^i^ and charge current N ^ from their local equilibrium forms. For example, in the local rest frame 
the stress tensor reads at leading order in gradients: Tij = Peq{e)Sij — r]{diUj + diUj + ^SijdkUkUj) — C^ij^ ' u. 
In this section, we will briefly review methods for calculating these transport coefficients (especially shear 
and bulk viscosity) in both weak and strong coupling regimes. 

For a weakly coupled system, the transport coefficients can be calculated from kinetic theory or from 
linear response theory. The kinetic theory approach [154ul^5lll56lll571IT58, 159 starts from local equilibrium 
distribution functions and expands the full distribution function order by order in terms of gradients of the 
four-velocity, temperature and chemical potential. The viscosity coefficients are associated with the first 
order terms in velocity gradients, and can be determined from the Boltzmann equations if the collision term 
is known explicitly (they depend on the transport cross section in the collision term) . 

In linear response theory, the transport coefficients can be rigorously expressed by Kubo formulae [160j , 
that relate the transport coefficients to the slope of associated spectral functions at zero frequency. For 
example, 

, = lim hm C = i lim lim (^^^^^ , (1.8) 

where the spectral functions are given by the imaginary part of the retarded green functions, p{oJ, k) = 
2ImG'ij(w, k), for certain components of the energy-momentum tensor: 

G^i^-'''''{u,'k)^~i [ d'*a;e*(^*-'''')([T^''(x,i)T'"^(0,0)]). 



The Kubo formulae can be applied to both weakly and strongly coupled systems (see Chap. 1.5) and 
can be generalized for lattice simulations. In lattice QCD, the spectral function p{uj) is obtained from the 
Matsubara (imaginary time) Green function Ge, instead of the retarded Green function Gr. (Both of 
them have spectral representation in terms of p ) The calculation of shear and bulk viscosity on the 
lattice |1611ll62lll631ll64] starts with the calculation of the corresponding imaginary time correlators of the 

^^Gu is obtained from Ge by replacing iuJn with w ~ po + is: 

duj p(a;) ^ / \ ■ f '^'^ P{^) 



GE(i'^' 



n) = / : — , Gfl(po) = -M ; 



2tt u — po — ie 

Here uj„ = 2nnT(n £ Z) are bosonic Matsubara frequencies, and T is the temperature of the medium. 
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energy momentum tensor in Euclidean time < r < l/T, which is related to the spectral function by 

duj p{uj) _ [ duj cosh[ijj{T ~1/{2T))] 



27r sinh[cj/(2r)] 



V(a;). (1.9) 



In principle, one can obtain the spectral function piu) by inverting ea. (|1.9p . and then use eq. (|1.8p to cal- 
culate the transport coefficients. In practice, one can only compute a finite number of points for Ge{t) 
on the lattice, which makes the unique extraction of an analytical spectral function p{u}) impossible. Gen- 
erally, one makes an ansatz for /o(w) with a small number of parameters, motivated by (usually somewhat 
model-dependent) considerations on the expected behavior of p{lo) at small and large w, and then fits these 
parameters to the Monte Carlo data [renfTO^llTMITM] . 



1.4 Shear and bulk viscosity for QCD matter 
Shear viscosity 

Attempts to calculate the QGP shear viscosity using weakly coupled QCD started more than 20 years 
ago [16511166] . Using kinetic theory in the relaxation time approximation and using a simple perturbative 
estimate for the latter, one found that the shear viscosity of the QCD matter behaves as ~ aHn{i/a ) l^^^^ 
1166] . A first calculation of the leading logarithmic contribution from the Boltzmann equation was performed 
in [167j . with a complete result finally published in 156]. A full leading order calculation that also computed 
the coefficient under logarithm was performed by Arnold, Moore and Yaffe using effective kinetic theory in 
the hard thermal loop approximation ^157j . For a weakly coupled QGP with three massless quark flavors, 
the shear viscosity to entropy ratio is [157] : 

■n _ 5.12 

This result is shown as the solid line in the left panel of Fig. ll.lOi using the two-loop renormalization 
group expression for the running coupling g{T) |145j . (Note that using Ba ll. 101 with a running coupling is 
phenomenological and beyond the order of accuracy at which Ba ll. 101 was derived.) 

In principle, the QGP shear viscosity can also be non-perturbatively calculated by lattice QCD using 
Kubo formulae. However, such calculations are highly non-trivial in the standard Monte-Carlo simulations 
due to the large noise to signal ratio for the relevant operators and the ill-posed inversion problem for the 
spectral function p(w), given the finite number of data points for the Euclidean correlators |163j . A first 
attempt to calculate the shear viscosity of SU(3) gluonic matter on an S'^ x 4 lattice comes from the pioneering 
work of Karsch and Wyld [161!, using a three-parameter ansatz for the spectral function /o(w). This method 
was later implemented by Nakamura and Sakai to calculate both shear and bulk viscosity on larger lattices 
(16^^ X 8 and 24^ x 8) with improved Iwasaki action [162] . They found that the shear viscosity to entropy 
density ratio is less than one, 7]/ s < 1, and that the bulk viscosity is zero within errors for the temperature 
range 1.4 < T/Tc < 1.8. A new evaluation of the shear viscosity for SU(3) gluon dynamics was performed 
by Meyer |163j using a two-level algorithm [168j , which dramatically improves the statistical accuracy of the 
relevant Euclidean correlators. He derived a robust upper bound rj/s < 1.0 for the shear viscosity entropy 
ratio and estimated that: 

r 0.134(33) (r=i.65r,), 

\ 0.102(56) (r=1.24T,). ^ ' 

Early calculations of the shear viscosity for the hadronic matter started from a theory of massless pious 
in the low-energy chiral limit, giving [170] 
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Figure 1.10: left: 77/s for low temperature hadronic phase and high temperature QGP phase |145| . right: 
7]/ s for hadron resonance gas in chemical and kinetic equilibrium [169] . 

where /tt — 93 MeV is the pion decay constant. The ratio r//s diverges at zero temperature; together with 
the logarithmic perturbative results discussed above, 77/3 is seen to reach a minimum near the QGP phase 
transition. A more detailed calculation of 77/5 for hadronic matter with both pions and kaons can be found in 
Ref. |170j ; it goes beyond the chiral approximation and includes intermediate resonances such as the p meson 
(see also |171j ). The left panel of Fig. ll.lOl shows these two results, which qualitatively agree with each other 
at low temperature but gradually deviate from each other above T > 100 MeV due to kaon excitations. 

Recently, rj/ s for a hadron resonance gas including a large set of hadron species up to 2 GeV was extracted 
from UrQMD^n simulations, using the Kubo formula |169j . Fig. II. IOI shows the extracted rj/s as a function 
of temperature T for a chemically equilibrated hadron resonance gas with zero chemical potential. Below 
T < lOOMeV, rj/s quickly rises with decreasing temperature, in qualitative agreement with the chiral pion 
result jl70]. Between 100 MeV and 160 MeV, rj/s saturates at a value ^ 1. For a hadron resonance gas 
out of chemical equilibrium, a similar tendency was demonstrated in Fig. 4 of Ref. |169j . where it is shown 
that at fixed temperature T, a non-zero baryon chemical potential reduces rj/s. For pb/T = 3.0 — 3.4 
(an unrealistically large chemical potential for RHIC energies), r]/s can decrease to 0.3-0.4 near the phase 
transition, and the plateau structure between 100 MeV and 160 MeV is replaced by a tendency that rj/s 
continues to decrease with increasing T. 

It is worthwhile to point out that Ref. [172j argued that including Hagedorn stated can significantly 
reduce ry/s in the hadronic phase near Tc since the highly degenerate Hagedorn states could dramatically 
increase the entropy density. The authors of Ref. [172] investigated Hagedorn state effects on 77/3 for a gas of 
pions and nucleons [173j and for a hadron-resonance gas with excluded volume corrections [174j and found 
that near Tc, rj/s can be significantly reduced to values close to the KSS bound rj/s = 1/Att (see Chap. 1.5). 

Bulk viscosity 

At sufficiently high temperature, the equation of state of a massless QGP satisfies p = e/3 which, on the 
classical level, leads to a vanishing bulk viscosity for the weakly coupled QGP. However, quantum corrections 
break the conformal symmetry and result in a non-zero bulk viscosity. Detailed calculations from Arnold, 

UrQMD: Ultra-relativistic Quantum Molecular Dynamics, a hadron cascade model JJ^j^TjUiT^ ■ 

Hagedorn states are highly unstable, massive hadronic resonances that exist close to Tc- These states are not included in 
the above UrQMD simulation. 
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Figure 1.11: left: C/s for high temperature QGP and for massless (left) and massive (right) pions [153j . 



Dogan and Moore [175] showed: 



c 



AalT^ 
ln[^*/mi3] 



(1.13) 



Here, A, fi* and mjj depend on the number of quark flavors in QCD. For QCD with three massless quark 
flavors, A = 0.657, /x* = 7.77r and to|, = l.Sg^T^ . Relating Eq. (|1.13p with the shear viscosity calculated 
within weakly coupled QCD |156lll57j one finds an approximate, but simple relationship between shear and 
bulk viscosity: 



C ^ 157;(c^ - 1/3)2 



(1.14) 



While this relation does not hold exactly [175j . it gives the right order of magnitude of C/77. 

The bulk viscosity of interacting massless pions was calculated by Chen and Wang in the framework of 
kinetic theory, combined with chiral perturbation theory. They found that C/s monotonically increases with 
temperature below T = 120 MeV; - ?t- The bulk viscosity for massive pions was calculated by Prakash et 
al. in the 1990's, using kinetic theory with experimental elastic scattering cross sections as inputs [170til71) . 
In contrast to the case of massless pions, one finds that C/s is a decreasing function of temperature. The 
bulk viscosity to entropy density ratio | for massless and massive poins, together with the result for a weakly 
coupled QGP, are illustrated in Fig. 11.111 

The behavior of bulk viscosity near the QCD phase transition was first investigated by Paech and Pratt 
within the linear sigma model |176j . Their work showed qualitatively that C/s reaches a maximum near the 
phase transition. This behavior was confirmed in later research by Kharzeev et al. |177pi77j and Meyer [164j . 
respectively. Using low energy theorems, Kharzeev et al. connected the bulk viscosity with the lattice inter- 
action measure e — 3p and extracted a temperature dependent (/s from lattice data for pure SU(3) gluon 
dynamics [177j and for full QCD |178j . Meyer directly extracted (/s from lattice QCD simulations by calcu- 
lating the corresponding correlation function for a pure SU(3)gluon plasma. He showed that the peak value 
of C/s near Tc can reach as high as 0.73 ^164j (a value that is around 10 times larger than the string theory 
estimates from holographical models [179pi80) ). Both of these methods involve an ansatz for the spectral 
function. However, the parameterized spectral functions used by these authors were challenged by Moore 
and Saremi |181j . who examined the behavior of the spectral function with analytical methods both near 
the QCD phase transition and in the weak coupling regime. 
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1.5 Shear viscosity and bulk viscosity from iV = 4 SYM and the AdS/CFT 
correspondence 



The transport coefRcients calculated from weakly coupled QCD are valid at sufRciently high temperature 
where the running coupling constant becomes small. However, the temperature reached in relativistic heavy 
ion collisions is not very high (top RHIC energies: T ~ 350 MeV, and LHC energies: T ~ 600 MeV). 
It is thus questionable to directly apply the weakly coupled QCD results at LHC/RHIC temperatures. 
Some insights for this problem can be gained from studies of iV = 4 supersymmctric Yang-Mills theory 
(SYM), where the shear viscosity can be calculated in both the strong coupling ^182 . 183J and weak coupling 
regimes |184j . Without too much effort, one can parameterize the behavior at intermediate coupling by 
interpolating between the two limits |184j . 

For = 4 SYM theory with infinite Nc and large, but finite t'Hooft coupling g^Nc (strongly coupling 
regime), rj/s was calculated by Buchel, Liu and Starinets (using the gauge/gravity correspondence, see 
below). They found [183] 

rj 1 
s 47r 



^ , 135C(3) 



(1.15) 



8(2(72iV^)3/2 ■ 

The rj/ s N = A SYM theory in the weakly coupled regime was investigated by the McGill group using the 
kinetic theory approach. They found a much smaller 7]/ s than what was previously found in weakly coupled 
QCD at the same value of the coupling constant. However, good agreement can be achieved after re-scaling 
each result by a combination of Debye screening mass and the number of degrees of freedom in the theory. 
Extending this mapping between weakly coupled QCD and weakly coupled N=4 SYM to lower temperature 
or into the strongly coupled regime, one finds that as = 0.5 in QCD corresponds to A = 4.7 in N=4 SYM, 
where rj/s from N=4 SYM is still relatively large, of the order of 1 |184j . 

The finite temperature version of the gauge/gravity (AdS/CFT) correspondence jl85yi86) . which maps 
finite temperature gauge field theory at strong coupling onto weakly coupled gravity in a curved space with 
a black whole, provides a new method to study strongly coupled systems and helps us to gain insights for 
the properties of a strongly coupled qgfEI. One example is the shear viscosity in the strong coupling limit. 
In the language of black holes gravity, the shear viscosity, defined by the Kubo formula, is connected with 
the absorption cross section of the black hole: r) oc limi^^o ctbh = where the last equality comes from the 
general theorem on black holes |187j which states that the cross section (Jbh = A is equal to the horizon 
area A for a broad class of black holes. On the other hand, the horizon area A also represents the entropy 
of the black hole: Sbh = j^- This shows that the shear viscosity to entropy ratio -q/s is a constant |188j . 
The constant was found by Kovtun, Son, and Starinets [146) . who showed that 

(1.16) 

s 47r 

is a universal value for a large class of black holes or branes. Given that the corresponding dual gauge field 
theories are very different, they conjectured that l/47r is an absolute low bound for rj/s, now known as the 
KSS bouncIEI. 

The bulk viscosity vanishes in conformal field theories such as A'' = 4 SYM. To apply string theory 
techniques to the calculation of the bulk viscosity in the strongly coupled regime thus requires investigating 
non-conformal deformations of the original AdS/CFT correspondence. Some recent developments can be 
foimd in Ref [179|,I191[I192|,I193] . Based on various holographical model computations, Buchel proposed a 
lower bound for the bulk viscosity to entropy ratio [194j : 

^->2{\-cl)'^. (1.17) 

S 6 S 

^'^The AdS/CFT correspondence can not be directly applied to QCD since QCD is not a Conformal Field Theory, and its 
gravity dual (if it exists) is not known. 

Some possibilities to violate the KSS bound are discussed in il89i\19Ujl . 
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(Note the different dependence on i — from the weakly coupled result (jl.l4p V Using a scalar potential 
tuned to reproduce the lattice equation of state, Gubser et al. calculated the bulk viscosity via gravity duals 
of non-conformal gauge theories [179[|180j . They found that C/s rises sharply near Tc, with C/s\t^ — 0.05. 



1.6 Notation and outline of this thesis 

Throughout thesis we adopt units in which h — c = kB — The metric tensor is always taken to be 
(7^" — diag(+l, —1, —1, —1). A'^'^ = g^'^ — u^u'^ is the projector onto the space transverse to the fluid velocity 
w^, Af^^u^ — 0. The partial derivative 9'' can then be decomposed as: 

where, — A^^'^d^ is the transverse component of the partial derivative 9^, and D = u^^d^ is the corre- 
sponding longitudinal component. In the fluid rest frame, D reduces to the time derivation and reduces 
to the spacial gradient. One therefore also uses the denotation / = Df. 

We also frequently use the symmetric, anti-symmetric and the ( ) brackets in Chap. 2, following Ref. |195| 
[T96l[T97l[T98] : 



Ai,B^) = i (^A^A^ + A^A^ - ^A-'^A^,) 

Using the above notations, the commonly used local fluid rest frame variables in dissipative viscous 
hydrodynamics are expressed in terms of the energy momentum tensor T'^" , charge current and entropy 
current 5'^ as following (see Chap. 2.1 and Ref. |195[fT96..197J for details): 



n = 




net density of charge; 


= 


AiiN" 


net flow of charge; 


e = 




energy density; 


P + n = 


-1/3 A^, T^"^ 


p: thermal pressure, 11 : 


Tt'"' EE 




shear stress tensor; 


W>' = 


u^.T^^A^l 


energy flow; 






heat flow; 


s = 




entropy density; 




Ai^S" 


entropy flux. 



In the 1st and 2nd order formalisms for viscous hydrodynamics (Chap. 2) one also frequently encounters 
the following scalar and tensors constructed from the gradients of the four velocity : 

9 = d ■ u expansion rate, 

V<''m''> = cr^'' = i(V^u'' -I- Vu^) - iV^'^a^u" velocity stress tensor, 

fj^" = -VI^m'^I vorticity tensor. 

Note that the above notations are for Cartesian coordinates — {t,x,y,z). In Chap. 3, we change to 
curvilinear coordinates = (r, x, y, 77) for the convenience of describing a longitudinally boost invariant 
system (i.e. viscous hydrodynamics in 2+1 dimensions). To express the above variables in [T,x,y^rj) co- 
ordinates, one can notationally replace by m,n, but must also replace the partial derivative by the 
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covariant derivative dm, see Chap. 3.1 and Appendix A. 1 for detail. 
Outlines: 

This thesis focuses on viscous hydrodynamics for relativistic heavy ion coUisions. Chapter 2 outUnes 
the 2nd order formalism for viscous hydrodynamics, obtained from different approaches. Chapter 3 sets 
us up for numerical calculations in (2+l)-dimension, assuming exact longitudinal boost invariance. The 
discussion there includes explicit transport equations in (2+l)-dimension, the initial conditions, the EOS, 
final conditions and the transport coefficients used in viscous hydrodynamic simulations. Generic shear and 
bulk viscosity effects are studied in Chapter 4 and Chapter 6, respectively, by comparing runs with ideal 
and viscous hydrodynamics using identical initial and final conditions. Chapter 5 focuses on the system 
size dependence of the shear viscous effects and investigates the multiplicity scaling of the elliptic flow 
coefficient V2 for both ideal and viscous hydrodynamics. In Chapter 7, we report on some of the most recent 
developments in viscous hydrodynamics, including the search for the optimized I-S equations for numerical 
implementations and recent efforts on code verifications among different groups. Chapter 8 assesses the 
current uncertainties for extracting the QGP viscosity from experimental data and briefly comments on 
some future directions for viscous hydrodynamics. Chapter 9 briefly outlines other methods for estimating 
the QGP shear viscosity. Our conclusions are summarized in Chapter 10. 
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Chapter 2: Relativistic Viscous Hydrodynamics 

— the Formahsm 



In this chapter we introduce the formahsm for viscous (dissipative) hydrodynamics. Initial attempts to 
formulate relativistic dissipative fluid dynamics started as a relativistic generalization of the Navier-Stokes 
(N-S) equations [199[|200| . Unfortunately, it turned out that the relativistic N-S equations are unsuited for 
numerical implementation since they developed instabilities due to exponentially growing modes, in partic- 
ular at high frequencies, and violate causality (the N-S formalism will be briefly described in Chap. 2.2.). 
These difficulties are largely avoided in the 2nd order formalism developed 30 years ago by Israel and Stew- 
art [201lll98ll202j . which goes beyond the so-called 1st order Navier-Stokes approach by expanding the entropy 
current to 2nd order in dissipative flows, replacing the instantaneous identification of the dissipative flows 
with their driving forces multiplied by some transport coefficient (as is done in Navier-Stokes theory) by a ki- 
netic equation that evolves the dissipative flows rapidly but smoothly towards their Navier-Stokes limit. The 
deduction of the I-S equation from the entropy current expansion and the 2nd law of thermodynamics will be 
introduced in Chap. 2.3. The I-S formalism can also be derived from the Boltzmann equation [131.203,159j. 
by expanding the distribution function around local equilibrium. The most general form of 2nd order vis- 
cous hydrodynamics was derived in |204|,I159] for systems with conformal symmetry and in [203] for systems 
without such symmetry (i.e. which also feature bulk viscosity and heat conductivity). Some of this will be 
discussed in Chap. 2.5. All versions of the I-S formalism directly solve evolution equations for the dissipa- 
tive flows. In contrast, another type of 2nd-order formalism, developed by the Ottinger and Grmela (0-G 
formalism) |2051l206|l207j . uses evolution equations for auxiliary fields, rather than the dissipative flows (see 
Chap. 2.6). The standard dissipative flows can be approximately reconstructed from the auxiliary fields [208j . 



2.1 From ideal to viscous hydrodynamics 

Hydrodynamics is a macroscopic tool to describe the expansion of the QGP and hadronic matter. It 
starts from the conservation laws for the conserved charge currents [x) and the energy momentum tensor 
T'^'^(x) {x denotes the 4-dimensional space-time coordinates (i,x)) [209j : 

d^Ni^{x) = Q, i = l,...,fc; (2.1a) 
d^T^^''{x)^0. (2.1b) 

For simplicity, one sets k = \ and only considers the conserved net baryon number current. This leaves 
14 independent variables: 10 from the symmetric energy momentum tensor T^" and 4 from the charge 
current N^. However, this system of equations is unclosed since (|2.ip only offers 5 independent equations. 
To solve this problem, one needs either to reduce the number of independent variables through physically 
assumptions, or to provide more equations. 

Ideal hydrodynamics [26l[27] solves this problem via the first route. By assuming perfect local thermal 
equilibrium, one can decompose the charge current and energy momentum tensor as follows: 

TV^ = nu'' , (2.2a) 
T"" = eu^'u" - pA^", (A''"' = g"" - u'^u^) , (2.2b) 

where n.(a;), e(a;), p{x) are the local net baryon density, energy density and pressure, respectively, and uf^{x) 
is the 4-flow velocity which is time like and normalized: u^u^ = 1. The above ideal fluid decomposition 
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reduces the 14 independent variables into 6 (1 each for e, n and p and 3 independent components in u^). 
With one additional input - the equation of state (EOS) p = p{n, e) - the system of equations is closed 
and can be solved to simulate the evolution of the system. Using the fundamental thermodynamic identity 
e + p ^ Ts + /in, it is easy to show that, in the absence of shocks (i.e. discontinuity in e, n or p), Eqs. (j2.ip 
conserve entropy, S^S"^ = 0, with S'^ = su'^. 

In classical kinetic theory, and T^'' are associated with the microscopic phase-space distribution 
function f{x,p) as follows [20911210] : 

N^^Y.'^^j^P>'h{x,p), (2.3a) 

i 

i •' 

Here are the baryon charges per particle of species i (particles and anti-particles count as sperate species) . 
Plugging in the equilibrium distribution function fi^ eq{x,p) — ^[p.u{^)+^^\^) ] /T(^)^i (where p,i{x) and T{x) are 
the local chemical potential of particle species i and the temperature, respectively) into ea. (|2.3p . one directly 
obtains the ideal fluid decomposition (|2.2p . Local thermal equilibrium is thus the basic assumption behind 
ideal hydrodynamics. In kinetic theory, it requires that the microscopic mean free path is much smaller than 
the system size and the microscopic collision time scale is much shorter than the macroscopic evolution time 
scale. The concept of local thermal equilibrium is, however, more general and can be formulated without 
recourse to kinetic theory. Ideal fluid dynamics always applies if local thermal equilibrium is ensured. 

If microscopic processes are not fast enough to satisfy the above conditions, the system is no longer 
in local equilibrium. For a near-equilibrium system, the distribution function can be decomposed into an 
equilibrium part plus a small deviation: 

/(x,p) = Ug{x,p) + Sf{x,p) , {6f « /) (2.4) 

Putting eq. (|2.4p into eq. (|2.3p one finds 

N'' = + SN^", (2.5a) 
Tf"" = T^^'' + , (2.5b) 

where the dissipative flows SN^^ and JT^'' are generated by the non-equilibrium contribution Sf. By imposing 
the "Landau matching conditions" [200] for an arbitrary frame u^: Uf^6N^ = and u^ST^^'^u^ = 0, one 
associates N^^ and T^" with the equilibrium definitions of net baryon density n and energy density e and 
pressure p (n = N^u^, e = u^T^^'^Up). In other words, these matching conditions fix the temperature and 
chemical potential of the equilibrium distribution f^q in (|2.4p such that e and n are defined in terms of feq 
in the standard way. Then , T^'' can be fully decomposed as [209] : 

Nf" = nuf" + y^, (2.6a) 
T^"^ = eu'^u'' - pl^"" - HA'"' + TT^'' + W^'u" + Wu^", (2.6b) 

where y , 11, tt'*" and W^^ are called dissipative or viscous flows. More specifically, — A'^^iV^ describes a 
baryon flow in the local rest frame, and = A^Ti^^u" = !^-i-£-V + q^^ (where q^^ is the "heat flow vector") 
describes an energy flow in the local rest frame. , and q^ are all transverse to the frame u^: W^Ufi = 0, 
V^Ufi, = and g'^w^ = 0, so each of them has 3 independent components. Since q^ is defined through 
and ^ the 3 vectors leave 6 independent components altogether. 11 ~ — ■|A^^T'^'^ — p is the viscous 
bulk pressure, which contributes to the trace of energy momentum tensor, tt^'' = y{Mi'> ^ [^(A^'^A'''^ -\- 
A''°'A^'^) — iA^" A'^'^jTi-cr (where the expression {p,v) is the shorthand for traceless and transverse to 
and u„ as defined by the projector in square brackets) is called the shear stress tensor, -k^" is traceless 
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and transverse to m'^, tt^ = and tt'^'^u^ = 0, and it is symmetric (like T'^"); it thus has 5 independent 
components. 

The Landau matching conditions, however, leave the choice of frame unconstrained. Two commonly 
used frames are the Eckart frame |199| and the Landau frame |200j . The Eckart frame sets u'^ parallel to the 
charge flow A^'^. As a consequence, disappears and the energy flow vector reduces to the heat flow vector 
= . However, for a system with very small or vanishing net baryon number, which is approximately 
realized in nuclear collisions at top RHIC energy, this frame is ill defined |211j . Therefore, this thesis and 
other theoretical viscous hydrodynamics frameworks aimed at RHIC physics |195II196II197II124II212II213| adopt 
the second choice ~ the Landau frame. In the Landau frame, is taken parallel to the 4-velocity of energy 
flow (cx T^^'^Ujy), which leads to a zero value for W^. Then the tensor decomposition shown in eq. I|2.6p 
simplifies as follows: 

77 

Nf" = nu" qf", (2.7a) 

e + p 

T'"' ^eu'^u" - {p + n)A'"' + n'"'. (2.7b) 

For the full evolution of all components of T'^'' and N^, we need in addition to the 5 conservation law 
equations (|2.ip additional equations for q^, H and w'^'^ . 



2.2 Navier-Stokes formalism 

The Navier-Stokes (N-S) formalism comes from the relativistic generalization of the non-relativistic N-S 
equations, which impose linear relationships between the dissipative flows and the corresponding thermody- 
namic forces [I99 1 I200 ] : 



A 



-V^ (^) ^ AX^ (2.8a) 

Il = -Ce = cx, (2.8b) 
tt'''' = 2ryV<^u''> = 2r;X^^ (2.8c) 

The (positive) "transport" coefficients C, A and 77 are the bulk viscosity, heat conductivity and shear viscosity, 
respectively. The scalar, vector and tensor thermodynamical forces on the right hand side are explicitly given 
hyX = -e = -d-u, X" = ^-ii" = -^V'f andX'"^ = V^'^u"') = a^"" = i(V^u'^-|-V''M^)- iA'^^^S^w". 

Unfortunately the N-S equations violate causality [2141l215ll216j and lead to infinite speed of signal 
propagation: if the thermodynamic forces turn off suddenly, the dissipative fiows will also disappear in- 
stantaneously. This confiicts with the fact that a macroscopic process caused by microscopic scattering is 
delayed by a relaxation time comparable to the kinetic scattering time scale. Connected with the acausality 
problem in the N-S formalism are numerical instabilities 214,215,216 , which render it useless for numerical 
simulations. Both of these two problems are avoided in the 2nd order formalism which is introduced in the 
following sections. 



2.3 Israel-Stewart formalism from the 2nd law of thermodynamics 

The 2nd order formalism for relativistic dissipative fiuids was first obtained by Israel and Stewart in 
the late 1970's, generalizing earlier non-relativistic work by I. Miiller [217j . They gave two derivations, a 
macroscopic one based on the 2nd law of thermodynamics, derived by expanding the entropy current up 
to the 2nd order in deviations from local equilibrium, and a microscopic one based on a near-equilibrium 
expansion of the Boltzmannn equation. Detailed presentations can be found in the original articles of Israel 
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and Stewart pOTll^U^] and in recent articles by Muronga [ro5inM[TF71l212l[^ . who first brought this work 
to the attention of the RHIC community. 

In a relativistic equiUbrium system, the entropy current is written as: 



(2.9) 



where a and (3 are related to the local temperature T and chemical potential /ibya=y,/3=y. For a 
near-equilibrium system, one can generalize the entropy current by the following off-equilibrium expansion: 

S^" = p(a, |3)|3^' - a{N^g + SN") + p^{T^^ + 5T^"') + Q^'{5N^', .JT^"^), (2.10) 

where, in addition to the first order corrections 5N^ and 5T^'^ appearing in the middle terms, the last term 
Q^^ includes second and higher order corrections in terms of 5N^^ and ST'^" . After some algebra, the entropy 
production rate can be written as: 



(2.11) 



After rewriting (5A^'' and JT'"' in terms of the corresponding dissipative flows 11, q'' and tt^", expressing 9^ta 
and d^Pi/ through the thermodynamics forces X^^ and X^^'^ defined in eq. (|2.8p . this can be further recast 
into: 



Td.S^ = nX - q^'X^ + TT^-^X., + T^M^' > 0. 



(2.12) 



The inequality on the right implements the 2nd law of thermodynamics d ■ S > 0. Note that the first three 
terms on the r.h.s. are first order, while the last term is higher order in the dissipative flows. 

The flrst order theory |,199,.200J neglects the second order and higher order contributions and sets = 0. 
By postulating linear relationships between the dissipative flows and the thermodynamic forces with non- 
negative coefficients described by the N-S equations (|2.8p . the inequality 9 • S* > is automatically satisfied: 



n2 



2AT 



2-q 



> 0. 



(2.13) 



Note that, since u^q^ = 0, q^ is space-like, g^g^ < 0. The N-S equations thus corresponds to a 1st order 
expansion of the entropy current, which is why the Navier-Stokes formalism is known as 1st order viscous 
hydrodynamics. 

The 2nd order theory of Israel and Stewart |20ip202(ll95pi96j keeps , more precisely, it keeps all terms 
that are second order in the dissipative flows: 



-{m^~l3iq^qf, + l32Vn 



2T 



a^Uq^ 
T 



T ' 



(2.14) 



where /?o, P2 and ao, cei are phenomenological expansion coefficients. After a bit of algebra, the entropy 
production rate Td ■ S is now written as il96til97j 



n 



1 



+ l3aU+ -Td. ^u'^ n - aaV^q' 



T 



InT - - f3iqi_, 



1 



-^2 A 1 , Y7 



(2.15) 
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Again, one can ensure the second law of thermodynamics by writing Tdf^^S^ in the form of eq. (|2.13p . This 
leads to the following 2nd order viscous equations for the dissipative flows: 



n = - 



1 



n + C0 - ZngV^q^ + ncT9^(^) 



(2.16a) 
(2.16b) 
(2.16c) 



where /3o, Pi and /32 have been replaced by the relaxation times rn = CPo, Tq = XT Pi and r^r = 2?7/32, 
respectively, and the mixing coefficients ag, ai have also been rewritten into luq = Cq^Oj ^qii — XTao, 
IqTT = XTai and l^rq = 2?7ai. The relaxation times replace the instantaneous identification between dissipative 
flows and thermodynamics forces in the N-S equations. They ensure the causality of the theory and the 
numerical stability [21411218] as long as they are not too small [2181. The viscous hydrodynamics based on 
the I-S formalism is thus also known as causal viscous hydrodunamicl^ 

In this phenomenological approach, both the viscous coefficients and the corresponding relaxation times 
are arbitrary free parameters, which must be determined from other theoretical consideration or extracted 
from experimental data. 



2.4 Israel-Stewart formalism from kinetic theory 

The kinetic theory approach starts from a Taylor expansion of the distribution function / around its 
local equilibrium form: 

/ = /o [1 + (1 ± fo)Sf] , Sf « 1, (2.17) 
with, 6 fix, t, p) = e(x, t) + eA(x, t)p^ + eA.(x, i)/p" + 0(/). 

Note that the scalar, vector and tensor coefficients e, e\, e\i, are functions of space-time, which, through 
kinetic theory, can be related to the dissipative flows. After putting eq. ([2.17P into the kinetic definition 
p.3p of T'^" and integrating out the momentum degrees of freedoms, one flnds that speciflc combinations 
of e, eA, correspond to bulk pressure 11, heat flow q\ and shear stress tensor ttxv For simplicity, we 
consider a fluid with only shear viscosity, neglecting bulk viscosity and heat conductivity, which simplifles 
the expansion of 5f as 

(5/(x,i,p)=eA.(x,i)pV, (2.18) 

with eai/ being traceless. Then, the shear stress tensor tt^" is directly related to e^'^ after integrating out 
the momentum degrees of freedom in eq. ([2.3p . For a massless Boltzmann gas /o = cxp (— p^u^/T), one 

finds [m] 

(2.19) 



2r2(e+p) ' 

with small quantum statistical corrections for fermions or bosons |137| . 

In the kinetic theory approach, the macroscopic conservation laws 9^iV^ = and d^T^^'^ — and the 2nd 
order I-S equations can be obtained from the Boltzmann equation p^d^f{x,t,p) — C{x) after integrating 

The projectors A*"^ on the r.h.s. of 12.16\ b. c) ensure that the transversality to the flow u^^ of and vr^^ is preserved 
during the time evolution. 
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out the momentum degrees of freedom, J duj = J ^^^^f^^ , with different combinations of momentum as 
weighting: 



dujp''df,f = J divC, (2.20) 

dcj p^'p'^d^J ^ J duj v'^C , (2.21) 

J dcj p^'p°'p^^df,f ^ J duj p°'p'^C . (2.22) 

Here C{x) is the coUision term whose precise form is related to the interaction matrix elements between 
particles. For a theory with conserved charges and with interactions that conserve energy and momentum: 
/dw C = and j duj p"C = 0. Eqs. ([MQllMll) then lead to ^^,N^' = and S^T^" = 0, respectively. After 
employing the relaxation time approximation for the collision term: 

C = -p^u''^-^^. (2.23) 

and integrating out J duj in eq. (|2.22p one finds the 2nd order I-S equations [131ll204j : 



2rj 



T^y^^'u"^ , (2.24) 



where ^^i, is the vorticity tensor, il^;^ = A^A(^(?[^Wq]. In contrast to the macroscopic approach of Chap. 2.3, 
the relaxation time r^r and shear viscosity rj here are no longer independent from each other. One finds 
Ttt = ^ for a massless Boltzmann gas, ''V = = It'^^^ ^'^^ ^ massive Boltzmann gas, and r^r ~ 1.024-^ 
for a massless Bose-Einstein gas |131[ . Actually, the shear viscosity r] can also be self-consistently calculated 
in the kinetic theory approach, after explicitly writing out the collision term. For purely gluonic matter, 
where interactions are dominated by 2 ^ 2 scattering processes with HTL cross sections, one finds: 77 = 

Based on consistent expansions to the 2nd order in Knudsen numbeij^, Betz, Henkel and Rischke |203j 
deduced the complete set of I-S equations from the Boltzmann equation for the bulk pressure 11, shear stress 
tensor tt^'' and heat flow vector q'' via Grad's 14-Momcntmii method [219j : 

n = HNs-Tnli 

+ An, (? • Va + An„ tt^'^ct^^ , (2.25a) 

- A,,a'^'^g, + A,nnV'^a + A,^7r'^"V,a , (2.25b) 
+ 2 T^q q^^'ii''^ + 2 l-^q V^^g''^ + 2 tt^^Q^^^ - 2 77 (Jz tt^'' 61 

- 2r^7r^<^CT''>^ -2A^,g<''V''>a + 2A^nnCT^'' . (2.25c) 

They found that the form of the I-S equations is independent of the choice of frame (Eckart frame vs. Lan- 
dau frame), whereas the coefficients rn, t,, r^, rng, T^n, t^tt etc. are frame dependent and are complicated 

^"^ Knudsen number is defined as microscopic mean free path over macroscopic scale of hydrodynamics, K = l-mfpl ^hydro- 
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functions of temperature and chemical potential. The first terms on r.h.s of each equation, corresponding to 
the Navier-Stokes first order approximation in Chap. 2. 2, are of first order in the Knudsen number. All other 
term on r.h.s are of second order in the Knudsen number. Compared with eq. (|2.16|) from the macroscopic 
approach, the Knudsen number expansion yields more terms for the 2nd order viscous equations. Some 
of these terms do not contribute to entropy production and thus can not manifest themselves in the 2nd 
law of thermodynamics; some others are suspected to correspond to higher order corrections in the entropy 
expansion approach |203| . 



2.5 Israel-Stewart formalism from conformal symmetry constraints 

For conformally invariant theories, the action S[(l),g^'^], as a functional of the external metric g^'^, is 
invariant under Weyl transformations: —^ e^'^^g^i,. As a consequence, the energy momentum tensor 
T^'^ is traceless and transforms as T'"' — + ^Gujrpni' ^^(jgj. Weyl rescaling in 4 dimensions. Correspondingly, 
the shear stress tensor also transforms as 

^f^<^ e^-^T:f^r (2.26) 

The above Weyl rescaling for the shear stress tensor tt'"^ gives a constraint for the possible form of the 
viscous equations for tt'^'^, if one assumes that the conformal invariance is preserved throughout the evolution 
of the system |204j . It is easy to prove that the 1st order N-S equation tt^"^ — 27]V^^m'^^ naturally satisfies this 
constraint, since the shear viscosity and velocity tensor transform as rj — > e^^rj and V^^u^'^ — > e'^'^V^'^u''^ 
respectively. 

The 2nd order viscous equations can be generalized from the 1st order one by adding 2nd order corrections: 
TT^" — 2?7V^^m'^^ + (2nd order corrections), satisfying the conformal constraint (|2.26p and the transversality 
and tracelessness properties: tt'^'^m^ — and tt^ = [204]. The last two constraints allow for 8 possible 2nd 
order corrections [204]: 

D^^lneD'^Mne, L'^^'D'^^ In e, V^^u-^^ (V^u") , P^^^ V^'^u'^ g^s^^^u^^i 

P^p'^^^u'^g^s^P^ P^^p n-^-f g-fSnP\ u^R''<^'•'^^us, R^'"'^ (2.27) 

where R^t^'^l^ is the Riemann tensor and R^'^ is the Ricci tensor. However, only 5 combinations of the above 
terms transform homogeneously under Weyl transformation. This determines the form of the 2nd order 
viscous equations, leaving 5 independent coefficients |204j : 



AJJA^L>7r"'3 + ^7r^''(Vau") 



3 



.h^^(^^^^^)>^ _ ^vr<^17''>^ - ^r!<^r!'^>^ . (2.28) 
2?]^ 2?7 2 

Note that the k, term is related to curved spaces and vanishes in flat space. The conformal symmetry con- 
straint requires that A^A^Dtt"'' must combine with the term |7r^''(VQ,u") such that the additional deriva- 
tives of the fmiction lli{x) in the Weyl transformation cancel between the two terms. The term ^tt*^^ {V au") 
term can be written in different ways as discussed in Chap. 2.7 below. 

While the conformal symmetry constraint approach provides the general form of the 2nd order viscous 
equations, it does not give the 2nd order coefficients Ttt, k, Ai, A2 and A3. They must be determined through 
explicit calculation. Such calculations were performed for weakly coupled QCD by York and Moore [159] 



and for strongly coupled N=4 SYM theory in [204j. Table [2T] summarizes the results. It is reasonable to 
assume that the strong and weak coupling results bracket the value for these coefficients in QCD at realistic 
values for the coupling strength. 
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Approaches 


Ttz Ai 


A2 


A3 


K, 


1) N=4 SYM 


2~ln2 I) 
2ttT 2-kT 







V 
ttT 


2) Kinetic theory 


5.0...5.9>) 4.1...5.2jj^ 

T s T s 


-2r]T^ 








Table 2.1: 2nd order Israel-Stewart Coefficients in the strong coupling [204] and weakly coupling [159] limits. 



2.6 Ottinger-Grmela formalism 

The TS formalism discussed in Chap 2.3-2.5 directly deals with the dissipative flows (tt^'', 11 and q^) 
and their viscous equations. In contrast, the Ottinger-Grmela formalism starts with an auxiliary field c^'^ 
and its evolution equation (0-G equation). The details for the Ottinger-Grmela formalism can be found 
in [205ll206ll207j . In numerical implementation, one solves the 0-G equation first and then relates the 
auxiliary field c^u to the commonly used dissipative flows in the subsequent spectra calculations [137] . 

If taking the 0-G equation and directly replacing the auxiliary filed c^'^ by the dissipative flows, one 
obtains an effective TS equation in the Ottinger-Grmela approach. For vanishing bulk viscosity and heat 
conductivity, the effective I-S equations for shear stress tensor ■n'^" was deduced by Dusling and Teaney, and 
has the following form [^[Mll^ : 



2?7V^^u 



6 6a 



— TT^'^ATr''^^ - r^7r<^A^^"'^\ (2.29) 

■n 

Where, the third order and higher order terms have been dropped, a and (3 are free parameters here. Gen- 
erally, one uses a to fix the relaxation time r^r = rj/pa and sets /3 = a/3 to ensure conformal invariance (see 
eq. (p:^ ). 



2.7 Final comments on the state of the formalism 



Comparing the 2nd order equations (|2.16cp . (H^St), 1^?^ and (P?^ superficially, one finds as 

the common terms that are obtained in all of the different approaches only the following: 



a/3 



(2.30) 



This so-called simplified I-S equation^ was used in our early calculations |1251ll26j and in Chaudhuri's 
work |142lll43lll44j . However this equation does not preserve the conformal symmetry for a conformal 
fluid (i.e. a fluid that consists of massless degrees of freedom). To fix this problem one needs to add an 



additional term |T^7r'^''(VaM") as pointed out in Ref [204] : 



a/3 



:r.7r^''(V„w") 



(2.31) 



In Ref. [127] we called this equation the full I-S equation. The last term in eq. (|2.31[) is clearly illustrated in 
Eqs. (|2.24p . (|2.28p and (|2.29|i^4 At flrst sight, eq. (|2.16cp seems to miss this term, but this naive impression 
is incorrect. For a conformally symmetric fluid, the temperature T is the only scale in the problem and 



-Eg. I2.25\l also contains a term of 7r'"^6 = 7r^''(Vati"), but leaves an undetermined coefficient in front of it. 
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therefore rj ^ s ^ and ~ hence riT/T.^ ^ T'^. In this limit, The last term in eq. (|2.16cp can then 
be rewritten as: 

(^"') = - 0) = -l7r^''(V„u"). (2.32) 

The last equality follows from the 2nd order hydrodynamic equation [134] . showing that, up to terms of 
sub-leading order, the last terms in eq. (|2.16c[) and eq. (|2.3ip agree. 

The "full I-S" equation, with slight variations, such as using the replacement eq. p.32p . has been estab- 
lished as the standard viscous equation for numerical implementations |127[|130[|134[|135[|136[|141j . It has 
received additional support by the observation made in [127] that it appears to minimize the dependence of 
the hydrodynamic evolution on the value of the microscopic relaxation time r^. While the variations related 
to the use of eq. p.32p are equivalent for conformal systems, they differ in principle for systems with an 
EOS that breaks conformal invariance, e.g. through a phase transition. These differences turn out to be 
negligible in practice [127j . For a conformal theory with vanishing chemical potentials, 4 other terms can 
be added to the right hand side of the "full I-S" equation, with additional coefficients Ai, A2, A3, and k 
(k = in Minkowski space), as shown in eq. (|2.28p [204j . These terms include couplings to the vorticity 
tensor [204j which turns out to be small in heavy-ion collisions if the initial longitudinal velocity profile is 
boost invariant [134j . Even more terms arise in a kinetic theory derivation that does not assume conformal 
symmetry and includes the effects from bulk viscosity and heat conductivity, as shown by Eqs. (|2.16cl[2.25p . 
Their coefficients can be obtained from kinetic theory |203|,I159| at weak coupling or from the AdS/CFT 
correspondence at strong coupling [204] (see Table |2T|) . The so far accumulated numerical evidence |127pi35j 
(see Chap. 7.3 and 7.4) suggests that, except for the last terms in Eqs. (|2.3ip . all other second order terms 
are unimportant in practice, but a systematic study that confirms this beyond doubt remains outstanding. 



28 



Chapter 3: Causal Viscous Hydrodynamics 

in 2+1 Dimensions 



Although the Israel-Stewart 2nd order formalism for relativistic dissipative fluid dynamics was established 
thirty years ago, its numerical implementations for describing relativistic heavy ion collisions started only 
quite recently [195pi96|ll97|ll24ll212|,l213j . Currently all of the available viscous hydrodynamics codes are 
restricted to either 1+1 dimensions [13ll[l32l[l33] or 2+1 dimensions [T25l[TMl[T271[T3ni[TM[T85l[TMl[T87l 
I140[ll42lll43|,ll44] . Here the first number indicates the number of spatial dimensions in which the code 
solves numerically for the evolution of the hydrodynamic fields, and the "+1" stand for time. All existing 
codes assume longitudinal boost invariance, which allows to treat the longitudinal expansion along the beam 
direction analytically. This is a good approximation at RHIC and LHC energies for particle production near 
midrapidity. (l+l)-d codes assume additionally azimuthal symmetry around the beam axis, allowing spatial 
expansion in the radial direction. With them, we can only simulate heavy ion collisions with zero impact 
parameter. To describe non-central collisions, and in particular anisotropic (elliptic) flow, requires a (2+l)-d 
code, which allows for anisotropic expansion in the two dimensions transverse to the beam. 

During the past years, the OSU group developed a (2+l)-dimensional viscous hydrodynamic code to 
describe the space-time evolution of the QGP and subsequent hadronic matter in the two dimensional trans- 
verse plane - VISH2+1 (Viscous Israel-Stewart Hydrodynamics in 2+1 dimensions) jl25[|126[|127[|13d] . In 
this chapter, we will summarize the crucial ingredients for the numerical calculations. These include the 
explicit form of the viscous hydrodynamic equations in 2+1 dimensions (Chap. 3.1), the initial conditions 
(Chap. 3.2), the equation of state (Chap. 3.3) and the freeze-out condition, including the freeze-out proce- 
dure and the calculation of spectra (Chap. 3.4). We also discuss the inputs for shear viscosity, bulk viscosity 
and the relaxation times (Chap. 3.5). 

3.1 Viscous hydrodynamic equations in 2+1 dimensions 

For ease of numerical implementation, the current viscous hydrodynamic calculations focus on boost- 
invariant systems, realized by assuming a specific "scaling" velocity profile = z/t along the beam direction. 
As shown by Bjorken |222j . this profile is a solution of the hydrodynamic equations (ideal or viscous) if the 
initial conditions are independent of the longitudinal reference frame (boost invariance), i.e. do not depend on 
space-time rapidity r/. After implementing the Bjorken approximation, the (3+l)-d viscous hydrodynamics 
reduces to (2+l)-d viscous hydrodynamics with boost invariance. Currently, all of the existing (2+l)-d 
viscous hydrodynamic codes |125[|126[|127[|135[|136[|137j also assume zero net baryon density and zero heat 
conductivity. These conditions are approximately realized in experiments at top RHIC and LHC energies 
and simplify the implementation further by eliminating the need to solve for the flows of baryon number and 
heat. 

Longitudinally boost-invariant systems are conveniently described in curvilinear coordinates = (t, x, y, rj), 
where r = V t'^—z'^ is the longitudinal proper time, V = ^^'^{jz^) is the space-time rapidity, and {x,y) are 
the usual Cartesian coordinates in the plane transverse to the beam direction z. In this coordinate system, 
the transport equations for the full energy momentum tensor T^'^ are written as |1241I126] : 



drT^^ + d,{v.,Tn + dyivyT^n = 

drf^^ + d^{v.,f-^) + dyiVyT^^) = 

drf^y + d^iu^f^y) + dyivyf^y) = s^y. 



(3.1a) 
(3.1b) 
(3.1c) 
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Here T"" = T(ro™"+7r""-nA""), Tg™" = eu"M"-pA™" being the ideal fluid contribution, = (w^, m^, m^', 0) 
7j.(1, i'x, Vy,0) is the flow proflle (with 7^ — ^ ), and g™" = diag(l, — 1, — 1, — 1/t^) is the metric ten- 

sor for our coordinate system. The source terms 5™" on the right hand side of Eqs. (|3.ip are given exphcitly 
as 

S^^ = -{p + n) - t^tt"" - ra:r(pwx+nw^+7r^^-wa;7r^^) 

- ra„(pwy+nt>j,+7r^^^-t;y7r^^), (3.2a) 
= -Ta:,(p+n+7r^^-W:,7r^^) - Ta2,(7r^'^-Wj,7r^^), (3.2b) 

s'^y = -Td^iir'^y-v^Tr^y) - Tdyip+u+iryy-vyTT^y). (3.2c) 

The transport equations for the shear pressure tensor and bulk pressure in 2+1 dimensions are written 
as [I261IT27] : 

1 1 ??T / r 

/^mn o„~mnA f^.m~nk \ ^,n~mk\ t^^, -ran I i I 



^^™„ = (^f™"~2?7CT"") - (u"^f"''+u"r"'')i:'ufe - -if^'-i-dfe ^w*^ , (3.3a) 

2 VTyT / 

i^n = -— (n-C0)-in^4. f^^i'^V (3.3b) 



rn 2 rn 

Here, we have written out the explicit (2+l)-d form for the full "I-S equations" described in Chap. 2.7. 
The expressions for (t"*" and tt™" are found in Eqs. (jA.2| IA.3|) in Appendix lA. 2.11 they differ from tt"™ in 
Eqs. p.2all3.2c|) and tr™" given in Ref. |124| by a Jacobian factor in the (?777)-component: tt'''' =t'^tt^^, 
^vv — j-'ig-m^ This factor arises from the curved metric where the local time derivative D = u™'dm must be 
evaluated using covariant derivatives d,n- Since = 0, no such extra Jacobian terms arise in the derivative 
Duk in the flrst line of Eq. (|3.3p . In more detail, D — u^dr +u^dx +uydy, 9 = d ■ u = drU^ + dxu'^ + dyU^ + ^ 
and cr™" = V<'"u"> =i(V"u"+V"u'")-iA""6l (with V^^A^'d,)- 

Even though several components of the symmetric shear pressure tensor tt™" are redundant [124] on 
account of its tracelessness and transversality to the flow velocity u™, VISH2+1 propagates all 7 non-zero 
components and uses the tracelessness and transversality conditions as checks of the numerical accuracy [125j . 
We flnd them to be satisfied with an accuracy of better than 1 — 2% everywhere except for the fireball edge 
where the tt™" are very small and the error on the transversality and tracelessness constraints can become 
as large as 5%. 



3.2 Initial conditions 

The initialization of a hydrodynamic simulation requires a starting time tq and initial profiles for the 
energy momentum tensor T'^'^, T'^^ and T"^^, which are given by the initializations for the energy density, ve- 
locity and stress stress tensor tt™". Through all of this thesis, we set tq = 0.6 fm/c, following the "standard" 
thermalization time used in most ideal hydrodynamic simulations 26J, and we use zero initial transverse 
flow velocity. We have not explored the need for rescaling tq when including viscosity; this awaits a careful 
comparison with the experimental data. The initializations for energy density and tt™" profiles are described 
below. 



3.2.1 Initializations for the energy density 
Glauber model initialization 

A simple Glauber model initialization assumes that the initial energy density in the transverse plane is 
proportional to the wounded nucleon density |223pi26) : 
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eoix, y; h) = Kn^;^^{x, y; b) 



(3.4) 



K^TA{x+^,y) 
+ TB{x-^,y) 



1-1- 



1 - ctTb {x-^,y) 



1-1 



B 



l-aTA{x+l,r^^^ 



A 



Here a is the total inelastic nucleon-nucleon cross section for which we take cr = 40 mb 

-oo 



clear thickness function of the incoming nucleus A or B, defined as Ta{x, y) = / 

Po 



Ta.b is the nu- 
dzpA{x,y,z); pA{x,y,z) 

is the nuclear density given by a Woods-Saxon profile: pAir) — i_)_gxp [ (r-fl^)/^ ] • ^^^^ i?cu=4.2fm, 
^ = 0.596 fm for Cu+Cu collisions and Rau = 6.37 fm, ^ = 0.56 fm for Au+Au collisions (these two parameter 
sets for Au and Cu nuclei correspond to a nuclear density po = 0.17fm^'^). The proportionality constant 
K does not depend on collision centrality but on collision energy; it fixes the overall scale of the initial 
energy density and, via the associated entropy, the final hadron multiplicity to which it must be fitted as 
a function of collision energy. The energy density profile for specific collision energy is therefore normal- 
ized by a parameter eo(T = tq, r = = 0) giving the peak energy density in the center of the fireball for 
central collisions (impact parameter b = 0). For central Au-|-Au collisions at top RHIC energies, one sets 
Co = e(0, 0; 6=0) = 30 GeV/fm'^ and tq = 0.6 fm/c to reproduce the final multiplicity in ideal hydrodynamic 
calculations. In the viscous hydrodynamics comparison runs, we use the same initialization as for ideal 
hydrodynamics. Viscous entropy production then leads to slightly larger final multiplicities than in ideal 
hydrodynamics. Again, we leave a corresponding retuning of initial conditions to a careful comparison study 
with experimental data. The Glauber initialization p.4p does not correctly reproduce the measured central- 
ity dependence of dNch/dy |224j . To fix this problem, more sophisticated initialization schemes have been 
developed [224j . where one uses a superposition of wounded nucleon and binary collision densities. We also 
leave this for future data comparisons. 

Color glass condensate (CGC) initialization 



The CGC initialization based on the Kharzeev-Levin-Nardi (KLN) approach [22511226] and its more 
recent fKLN improvement |227l has been applied in earlier ideal hydrodynamic calculations j2281l229[l230| 
and recently also in viscous hydrodynamic calculations [13511136] . Details of the CGC initialization can be 
found in Ref. [MHUHj- 




Figure 3.1: The initial eccentricity e^, from Glauber and CGC initializations, for Au-|-Au collisions at different 
centrality [135] . 
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Compared with the Glauber initialization, the CGC initialization gives a more "plateau-like" initial pro- 
file for the energy density |230j . which leads to larger initial eccentricity Ex = (l^'i^yi-^ than for the Glauber 
model. Fig l3. II shows the initial eccentricity Sx from the Glauber and fKLN initializations |135| . One finds 
that the CGC initialization predicts 20 — 50% larger initial eccentricity than the Glauber one, especially at 
small impact parameters, where the eccentricity is small but the discrepancy between models is large [231) . 

3.2.2 Initializations for tt"*" and 11 

Lacking a microscopic dynamical theory for the early pre-equilibrium stage, initializing the viscous pres- 
sure tensor tt™" requires some guess-work. We here explore two options: (i) zero initialization, which 
sets ^^"" = at initial time tq [HIIIISIIMI ■ (ii)(Navier- Stokes) N-S initialization: tt^" = 277crJ^" , where 
the shear tensor cr^" is calculated from the initial velocity profile u" = (l,0,0,0) [l25l [^61 [T27l [T35l 
I136[|137| . The second option is the default choice for most of the results shown in this thesis. It gives 
T^TTg'' — = ^^ttq^ — i.e. a negative contribution to the longitudinal pressure and a positive con- 

tribution to the transverse pressure. 



3.3 The equation of state (EOS) 

The equation of state (EoS), which relates the pressure p to energy density e and net baryon density ub'- 
p — p{7iB, e) (or equivalently relates it to temperature T and net baryon chemical potential /i^, p ~ p{T, fj-s)) 
is a necessary input for both ideal and viscous hydrodynamic simulations. In this section, we will describe 
the different EoS used in this thesis [1261ll27j . Currently, all of the existing viscous hydrodynamics calcu- 
lation neglect net baryon density and set ns = 0. The EoS described below thus reduces to a simple one 
dimensional EOS: p = p{e). 

EOS I : 

EOS I models a non-interacting gas of massless quarks and gluons, with p — ^e. It has no phase transition. 
Where needed, the temperature is extracted from the energy density via the relation e = (l6 + ^Nf^ fo (SHp'' 
corresponding to a chemically equilibrated QGP with Nf — 2.5 effective massless quark flavors |223) . 

EOS Q : 

EOS Q 223J connects a noninteracting massless QGP gas to a chemically equilibrated hadron resonance 
gas through a first order phase transition. The EOS in the QGP phase is defined by p= ^e—^B (i.e. 
c'^ = ^ — i). The vacuum energy (bag constant) B^^'^ = 230 MeV is a parameter that is adjusted to yield a 
critical temperature Tc — 164 MeV. The hadron resonance gas below Tc can be approximately characterized 
by the relation p = 0.15e (i.e. =0.15) [223j . The two sides are matched through a Maxwell construction, 
yielding a relatively large latent heat Aeiat = 1.15 GeV/fm^. For energy densities between en = 0.45 GeV/fm^ 
and cq = 1.6 GeV/im.^ one has a mixed phase with constant pressure (i.e. = 0). 

SM-EOS Q : 

SM-EOS Q is a smoothed version of EOS Q |126j . The discontinuous jumps of in EOS Q from a value of 
1/3 to at eg and back from to 0.15 at en generate propagating numerical errors in VISH2+1 which grow 
with time and cause problems. We avoid these by smoothing the function c^(e) with a Fermi distribution of 
width Se — 0.1 GeV/fm'^ centered at e — eq and another one of width 6e = 0.02 GeV/fm^ centered at e = cr. 
Both the original EOS Q and our smoothed version SM-EOS Q are shown in Figure 13.21 A comparison of 
simulations using ideal hydrodynamics with EOS Q and SM-EOS Q is given in Appendix IA.3.11 It gives an 
idea of the magnitude of smoothing effects on the ideal fluid evolution of elliptic flow. 
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Figure 3.2: The equation of state. 
Panel (a) shows the pressure p as a 
function of energy density e and (in 
the inset) the squared speed of sound 
= ^ as a function of temperature 
T, for EOS I, EOS Q, SM-EOS Q and 
EOS L. Panel (b) shows as a func- 
tion of energy density e for different 
EOS. 
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EOS L : 

EOS L [127] matches the hadron resonance gas below smoothly in a rapid cross-over transition to 
lattice QCD data [232j above Tc. For the fit, the lattice data were plotted in the form p{e), interpolated and 
then smoothly joined to the p{e) curve of the HRG with a cross-over transition near ~ 175MeV. As can 
be seen in the upper panel of Fig. 13.21 in the inset, our procedure is not fully thcrmodynamically consistent 
and leads to a somewhat different temperature dependence of below than for EOS Q and SM-EOS Q. 
Since this only affects the flow dynamics below our decoupling temperature of Tdec = 130 MeV, we have not 
put any effort into correcting thiJ^. 

PCE-EOS for the HRG phase: 

EOS Q assumes that chemical and thermal freeze-out happen at the same temperature Tth ~ 130 MeV [Ml 
I234j . However, the experimental data indicate that chemical freeze-out happens earlier near Tc ^ 165 MeV [Ml 
199] . The unrealistic implementation of chemical freeze-out in EOS Q leads to wrong predictions for the par- 
ticle ratios which need to be fixed by hand using normalization factors. 

Microscopically, chemical and thermal freeze-out are related to inelastic and total (elastic -I- quasi-elastic 
-I- inelastic) scattering rates between particles, respectively. The cross sections of the inelastic, particle 
number changing processes, are smaller than the elastic and quasi-elastic ones, which leads to an earlier 
chemical than thermal freeze-out |32j . A complete chemical freeze-out corresponds to a picture where 
all hadron numbers are fixed below Tch- However, quasi-elastic resonant scatterings constantly change 

The EOS L here is different from the cross-over EOS used in Ref which was constructed by connecting a hadron 

resonance gas EOS with the QGP EOS from a high-order weak-coupling QCD calculation J233f . 



33 



resonances and their daughter particles (e.g. tttt —> p tttt, ttN A ^ ttN, ttK K* ttK, etc.), such 
that the resonances are in relative chemical equilibrium with their decay products and their abundances 
only freeze out at thermal decoupling. This partial chemical equilibrium picture is naturally incorporated 
in any hadron cascade model, that solves the Boltzmann equations with elastic, quasi-elastic and inelastic 
cross sections for different hadron species |471I38[|37] . In pure hydrodynamic simulations it is realized by 
introducing an EOS describing partially chemical equilibrium (PCE-EOS) with effective chemical potentials 
adjusted to conserve the relative hadron abundances after resonance decays in the HRG phase |235[|236|I237[ 
1238] . Although the effective chemical potentials only slightly change the equation of state p{e), they affect 
significantly the relation between temperature and energy density T(e) in the HRG phase, since a larger 
portion of the energy density is stored in particle rest masses, reducing the thermal energy and temperature. 
As a result, the typical freeze-out temperature drops from 130 MeV (EOS Q) to 100 MeV (PCE-EOS) for 
the same freeze-out energy density edec = 0.085 GeV/fm^. 

Currently, the PCE-EOS has been well studied and implemented in ideal hydrodynamic simulations, but 
has not been applied to any viscous hydrodynamic simulations. However, the chemical content in the HRG 
phase at freeze-out will affect the extracted value of the QGP shear viscosity. This will be further discussed 
in Chap. 8. 

3.4 Freeze-out procedure and calculation of spectra 

The hadron spectra are computed from the hydrodynamic output via a modified Cooper-Frye procedure 
|239j . We here compute spectra only for directly emitted particles and do not include feeddown from 
resonance decays after freeze-out. We first determine the freeze-out surface S(a;), by postulating (as common 
in hydrodynamic studies) that freeze-out from a thermalized fluid to free-streaming, non-interacting particles 
happens suddenly when the temperature drops below a critical value. As in the ideal fluid case with EOS Q 
[223j we choose Tdcc = 130 MeV. The particle spectrum is then computed as an integral over this surface [239j . 



where gi is the spin-isospin degeneracy factor for particle species i, (Pijf^{x) is the outward-pointing surface 
normal vector on the decoupling surface Y,(x) at point cc, which for boost-invariant freeze-out at longitudinal 
proper time T/(r) reads 



(with r = (x,y) = [r cos (/), r sin (/)) denoting the transverse position vector), and fi{x,p) is the local distribu- 
tion function for particle species i, computed from the hydrodynamic output. Equation p.5p generalizes the 
usual Cooper-Frye prescription for ideal fluid dynamics 239 by accounting for the fact that in a viscous fluid 
the local distribution function is never exactly in local equilibrium, but deviates from its local equilibrium 
form by small terms proportional to the non-equilibrium viscous flows |240lll26lll31| . Both contributions can 
be extracted from hydrodynamic output along the freeze-out surface. The equilibrium contribution is 




(3.5) 



p- d^alx) = [rriT cosh{y—r]) p± ■ V^r/(r)] x Tf{r)rdr d(j)dri, 



(3.6) 




(3.7) 



where the exponent is computed from the temperature T{x) and hydrodynamic flow velocity 
it^ = 7x(cosh77, Wj;, Wj^, sinh?]) along the surface S(x): 



p ■ u{x) =7_L [mT cosh(2/-77) - pxVx - PyVy]. 



(3.8) 
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Here rriT = ^/p^+mj is the particle's transverse mass. 

As long as only shear viscosity is implemente d^. the viscous deviation from local equilibrium is given 

by [llollm] 

The viscous correction is proportional to tt^'^^x) on the freeze-out surface (normalized by the equilibrium 
enthalpy e+p) and increases quadratically with the particle's momentum (normalized by the temperature 
T). At large pT, the viscous correction can exceed the equilibrium contribution, indicating a breakdown of 
viscous hydrodynamics. In that domain, particle spectra can not be reliably computed with viscous fluid 
dynamics. The limit of applicability depends on the actual value of tt^'^ /{e+p) and thus on the specific 
dynamical conditions encountered in the heavy-ion collision. 

In Eq. (13. 9p . the final expression for viscous correction to the spectrum are written explicitly as: 

Pi^lP^t^^" = TO^(cosh^(?/-7])7r'^^ + sinh^(t/-?7)r^7r'''') 
- 2mT cosh(?/-?7) (pj,7r^^ -I- PyTr^^) 

+ {pI^^'-' + 2p^py^-y +plT:yy). (3.10) 



Due to longitudinal boost-invariance, the integration over space-time rapidity r] in Eq. (13. 5p can be done 
analytically, resulting in a series of contributions involving modified Bessel functions [32LI132] . VISH2-I-1 
does not exploit this possibility and instead performs this and all other integrations for the spectra numeri- 
cally [IMj . 

Once the spectrum p.Sp has been computed, a Fourier decomposition with respect to the azimuthal 
angle 0p yields the anisotropic flow coefficients. For collisions between equal spherical nuclei followed by 
longitudinally boost-invariant expansion of the collision fireball, only even-numbered coefficients contribute, 
the "elliptic flow" V2 being the largest and most important one: 

^^^'^ = dyprd;,d^} '^ = Y.^^^' + ^) ^'^-"^^^ + ■■■]■ (^-ll) 

In practice it is evaluated as the cos ( 2 0p) -moment of the final particle spectrum, 

f d4>r, cos(2fl!)„) f^-j-r 

v,{pr) ^ (cos(2<^p)) ^ ^ ' '^vUy,., ^^ ^3^^2) 

J dypTdpT d(j}p 



where, according to Eq. (j3.5[) . the particle spectrum is a sum of a local equilibrium and a non-equilibrium 
contribution (to be indicated symbolically as iV = N^q + SN) . 



3.5 Additional viscous inputs: shear viscosity, bulk viscosity and relaxation 
times 

The shear viscosity rj, bulk viscosity ( and their relaxation times and rn are free inputs in viscous 
hydrodynamic calculations. 

Although the future trend for viscous hydrodynamic calculations will be to input a temperature dependent 
77/s to consider the fact that 77/s is small in the QGP phase, reaches a minimum during the phase transition 
and rises again to much larger values in the hadronic phase (see Chap. 1.4 for details), all of the presently 
existing viscous hydrodynamic calculations |125pi26pi27pi35pi36pi37j (including those presented here) use a 

For the case with both shear and bulk viscosity, the expression for 5f can be found in Ref. \241^ - 
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Figure 3.3: Minimally constructed bulk viscosity to entropy ratio, C/s, as a function of temperature. 



constant r]/s as input: -q/s = C x j^, with C = 0, 1, 2, 3 ... is the minimal KSS bound from AdS/CFT 
[146j ). Typically, we set t^^ = except where mentioned otherwise. This expression is very close to 
the AdS/CFT prediction — ^"^^^ [242|,l243j . and is half of the kinetic theory prediction for massless 
Boltzmann particles [201|,ll31j . 

The bulk viscosity of the QCD matter is still under theoretical development (see Chap. 1.4, 1.5 for 
details). It is generally believed that C/s peaks near the phase transition [176|,I178| . However, one finds 
that the minimal peak value from AdS/CFT predictions |180pi79) is more than 10 times smaller than the 
one extracted from lattice QCD data [164lll77lll78j^ . Things become even more complicated for the bulk 
viscosity in the hadronic phase [153] . When neglecting hadron masses, chiral perturbation theory predicts 
a rising (C/s)(T) below Tc |244| . In contrast, a theory with non-zero hadron masses shows a decreasing 
{C/s){T) below Tc |1701ll71j . Considering these theoretical uncertainties, C,/s here is treated as a free input. 
We concentrate on the bulk viscosity effects near the phase transition, but totally neglect bulk viscosity 
effects in the hadronic matter. Fig. 13.31 shows the minimally constructed (C/'S)„j„, which is obtained by 
connecting the minimal strong coupling AdS/CFT result C/s = 2(?7/s)(l/3 — c^) [194] above Tc through 
a Gaussian function peaked at Tc with a zero value in the hadronic phase (the speed of sound, c^(T), is 
evaluated from the same lattice QCD data |232j that are used in our EOS L). To simulate effects from larger 
bulk viscosity, we multiply the entire function (C/s)min(?c) by a constant C > 1. For the bulk relaxation 
time rn, there are no standard theoretical results (see, however, the recent work [245] '). In Chap. 6.2, we 
set Tjj = = ^Jy just for the purpose of a qualitative comparison of shear and bulk viscous effects without 
claim of quantitative accuracy. In Chap. 6.3 and Chap. 6.4, we will investigate the effects of critical slowing 
down during the phase transition and corresponding bulk viscosity effects. There we will set Tn to be a) a 
temperature-independent constant that we change it from 0.1 fm/c to 10 fm/c; b) a temperature dependent 
function rn(T) = max[f • |(T'), 0.1 fm/c] with f = 120 fm/c. This choice implements phenomenologically the 
concept of critical slowing down; it yields Tn ~ 0.6 fm/c at T = 350 MeV and m « 5 fm/c at Tc. 



Refer to Ref. [ISl^ for a critical discussion of the lattice QCD approach. 
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Chapter 4: Generic Viscous Effects 

— Shear Viscosity 



4.1 Introduction 

In this chapter, we will study generic shear viscous effects by comparing ideal and viscous hydrodynamic 
runs. The results shown here are for Au+Au collisions and complement those from our two early viscous 
hydrodynamics papers |1251ll26j where we presented results for the smaller Cu+Cu systems. The choice of 
the Cu+Cu system was purely technical since the numerical grid in early VISH2+1 could not accommodate 
the larger Au+Au collision fireballs. This problem was solved in an up dated version of VISH2+1 with a 
new treatment of the boundary. However, since a similar detailed and systematic investigation for the shear 
viscous effects for Au+Au collisions have not been documented anywhere else, I decided to redo the graphs 
from [126| for Au+Au collisions and include them in this thesis, instead of using the old Cu+Cu graphs. 

As discussed in Ref. |127) . the shear viscous effects are qualitatively similar in Cu+Cu collisions and 
Au+Au collisions. However, on a quantitative level, the smaller Cu+Cu system shows larger viscous ef- 
fects |127j . This manifests itself through a much larger viscous V2 suppression for non-central Cu+Cu 
collisions than for non-central Au+Au collisions. The new Au+Au results were calculated with the updated 
VISH2+1, which solves the "full TS equation" rather than the simplified I-S equation used in our early 
papers |125[|126] . and organized along lines similar to Ref. |127j for comparison and later reference. As 
briefly mentioned in Chap. 2.7, the "full" I-S equation is preferred since it preserves conformal symmetry for 
a conformal fluid and reduces the dependence on the relaxation time r^r. Systematic numerical comparisons 
between the "full" and "simplified" I-S equations were presented in Ref. [l!27 , and will be discussed later in 
Chap. 7. 

Let me begin with a description of the default settings for the free parameters (see Chap. 3) used in the 
calculations for this chapter. 

For ideal and viscous hydrodynamic comparison runs, we use identical initial and final conditions, with 
the same standard parameters previously used in early ideal hydrodynamic simulations. Following Refs. |223| 
1234] . we use a Glauber initialization for the energy density and set eo(To,6 = 0,r = 0) = 30 GeV/fm'^ at 
To = 0.6 fm/c for central Au+Au collisions. The freeze-out surface is extracted using the AZHYDRO 
decoupling algorithm |246j with a constant decoupling temperature T^ec = 130 MeV, and the final particle 
spectra are calculated with the modified Cooper-Fye formula, given by eqs. ()3.5II3.10|) in Chap 3.4. In the 
viscous hydrodynamic calculations, we use the KSS value for the shear viscosity 77/s — I/Att |146j . and set 
bulk viscosity to zero: (/s = 0. For the relaxation time, we use r^r = irj/sT, and the shear stress tensor is 
initialized by the N-S initialization 7r'""'(To) = 2?7(t™", unless otherwise noted (e.g. in Chap. 4.4, where we 
explore the effects from varying the relaxation time and for shear stress tensor initialization) . 

To isolate effects introduced by the phase transition from generic shear viscous effects, we perform cal- 
culations with two different EOS, using EOS I for a pure massless quark-gluon gas without phase transition 
and SM-EOS Q for a more realistic EOS that includes a phase transition between QGP and hadron gas. 
Results from the lattice-QCD based EOS (EOS L), can be found in Chap. 6, where we also study the the 
bulk viscous effects. 
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Figure 4.1: Time evolution of the hydrodynamic source terms (|3.2aj|3.2c)) . averaged over the transverse 
plane, for central Au+Au collisions, calculated with EOS I in the left panel and with SM-EOS Q in the 
right panel. The smaller insets blow up the vertical scale to show more detail. The dashed blue lines are 
for ideal hydrodynamics with eo = 30GeV/fm^ and To = 0.6fm/c. Solid red lines show results from viscous 
hydrodynamics with identical initial conditions and 2 ^ ^ « 0.08, = ||; « 0.24 (200MoV^ ^^j^ rpj^g ^^^^^ 
in left panel indicate the freeze-out times at Tc = 130 MeV. The positive source terms drive the transverse 
expansion while the negative ones affect the longitudinal expansion. 



4.2 Hydrodynamic evolution 
Central collisions : 

Even without initial transverse flow, the N-S initialization tt™" = 2ria"^" generates three non-zero values 
for the components of the shear viscous pressure tensor: t'^tt^^ = tt^^ = 11^^ = J^, due to the boost- 
invariant longitudinal expansion. Inspection of the source terms in Eqs. p.2all3.2c| reveals that the initially 
negative t'^tt'^^ effectively reduces the longitudinal pressure, thus reducing the cooling rate, while the initially 
positive values of tt^^ and Tv^y effectively increase the transverse pressure and accelerate the development 
of transverse flow in x and y directions. As the fireball evolves, the shear viscous pressure tensor deviates 
from its N-S value (although for small relaxation time t^^ , the deviation remains small 0) , and the stress 
tensor cr™" receives additional contributions involving the transverse flow velocity and its derivatives (see 
Eq. (|A.3P ). which renders an analytic discussion of its effects on dynamics impractical. 

Figure 14.11 shows what one gets numerically. Plotted are the source terms p.2a|) and (j3.2bp , averaged 
over the transverse plane with the energy density as weight function, as a function of time, for evolution of 
central Au-|-Au collisions with two different equations of state, EOS I and SM-EOS Q. (In central collisions 
{\S^^\) = {\S'^y\).) One sees that the initially strong viscous reduction of the (negative) source term S'^'^ , 
which controls the cooling by longitudinal expansion, quickly disappears. This is due to a combination of 
effects: while the magnitude of t'^tt^^ decreases with time, its negative effects are further compensated by a 
growing positive contribution T{dx{pVx)+dy{pVy)) arising from the increasing transverse flow gradients. In 
contrast, the viscous increase of the (positive) transverse source term S"^^ persists much longer, until about 
6fm/c. After that time, however, the viscous correction switches sign (clearly visible in the upper inset 
in the right panel of Fig. 14.1b ) and turns negative, thus reducing the transverse acceleration at late times 
relative to the ideal fluid case. We can summarize these findings by stating that shear viscosity reduces 
longitudinal cooling mostly at early times while causing initially increased but later reduced acceleration 

These deviations from the N-S limit are investigated in Sec \7.3.7\ 
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Figure 4.2: Time evolution of the local temperature and local entropy density (down) in central Au+Au 
collisions, calculated with EOS I and SM-EOS Q (right), for the center of the fireball (r = 0, upper set of 
curves) and points at r = 3fm and r = 5fm (middle and lower set of curves). For clarity, the temperature 
curve for r = 3 fm and r = 5 fm have been lowered by constant offsets, as indicated. Same parameters as in 
Fig. 14.11 See text for discussion. 



in the transverse direction. Due to the general smallness of the viscous pressure tensor components at late 
times, the last-mentioned effect (reduced acceleration) is not very strong. 

The phase transition in SM-EOS Q is seen to cause an interesting non-monotonic behaviour of the time 
evolution of the source terms (right panel in Fig. 14. ip . leading to a transient increase of the viscous effects 
on the longitudinal source term while the system passes through the mixed phase. 

The viscous slowing of the cooling process at early times and the increased rate of cooling at later times 
due to accelerated transverse expansion are shown in Figure l4?2] bv the time evolution of local temperatures. 
For comparison we also show curves for boost-invariant longitudinal Bjorken expansion without transverse 
flow, labeled "(O-l-l)-d hydro". These are obtained with flat initial density profiles for the same value eg (no 
transverse gradients). The dotted green line in the left panel shows the well-known T^t^-^^^ behaviour of 
the Bjorken solution of ideal fluid dynamics [222j . modified in the right panel by the quark- hadron phase 
transition where the temperature stays constant in the mixed phase. 

The dash-dotted purple line shows the slower cooling in the viscous (0+l)-dimensional case |211j . due 
to reduced work done by the longitudinal pressure. The expansion is still boost-invariant a la Bjorken |222j 
(as it is for all other cases discussed in this thesis) , but viscous effects generate entropy, thereby keeping the 
temperature at all times higher than for the adiabatic case. The dashed blue (ideal) and solid red (viscous) 
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Figure 4.3: Surfaces of constant temperature T for central Au+Au collisions, evolved with EOS I (left panel) 
and SM-EOS Q (right panel). In each panel, results from viscous hydrodynamics (red lines) are directly 
compared with the corresponding ideal fluid evolution (blue lines). The right panel contains two isotherms 
for Tc = 164 MeV, one separating the mixed phase (MP) from the QGP at energy density eg = 1.6 GeV/fm^, 
the other separating it from the hadron resonance gas (HRG) at energy density ch = 0.45 GeV/fm'^. See text 
for discussion. 



lines for the azimuthally symmetric (l+l)-dimensional case show the additional cooling caused by transverse 
expansion. Again the cooling is initially slower in the viscous case (solid red), but at later times, due to faster 
build-up of transverse flow by the viscously increased transverse pressure, the viscous expansion is seen to 
cool the fireball center faster than ideal hydrodynamics. (Note also the drastic reduction of the lifetime of 
the mixed phase by transverse expansion; due to increased transverse flow and continued acceleration in the 
mixed phase from viscous pressure gradients, it is even more dramatic in the viscous than the ideal case.) 

The lower panels of Fig. 14.21 show the evolution of the entropy density. (In the QGP phase s ^ T^.) The 
double-logarithmic presentation emphasizes the effects of viscosity and transverse expansion on the power 
law s(r) ~ T~": One sees that the scaling of the ideal Bjorken solution is flattened by viscous effects, but 
steepened by transverse expansion. As is well-known, it takes a while (here about 6-8fm/c, depending on the 
EOS) until the transverse rarefaction wave reaches the fireball center and turns the initially 1-dimensional 
longitudinal expansion into a genuinely 3-dimensional one. When this happens, the power law s(r)^r~" 
changes from a = 1 in the ideal fluid case to a > 3 [26]. Here 3 is the dimensionality of space, and the fact that 
a becomes larger than 3 reflects relativistic Lorentz-contraction effects through the transverse-flow-related 
7_i_-factor that keeps increasing even at late times. In the viscous case, a changes from 1 to 3 sooner than 
for the ideal fluid, due to the faster growth of transverse flow. At very late times (not shown in Fig. ?? ) 
the s(r) curves for ideal and viscous hydrodynamics are almost perfectly parallel [126], indicating that very 
little entropy is produced during this late stage. 

In Fig. l4.3l we plot the evolution of temperature in t space, in the form of constant-T surfaces, obtained 
from ideal hydrodynamics (blue lines) and viscous hydrodynamics (red lines), respectively. Again the two 
panels compare the evolution with EOS I (left) to the one with SM-EOS Q (right). As already noted, at r = 
the viscous fluid cools initially more slowly (thereby giving somewhat longer life to the QGP phase) but later 
more rapidly (thereby freezing out earlier for SM-EOS Q). Comparing the two sets of temperature contours 
shown in the right panel of Fig. 14. 3i one sees that viscous effects tend to smoothen any structures related to 
the (first order) phase transition in SM-EOS Q. The reason for this is that, with the rapid change of the speed 
of sound at cither end of the mixed phase, the radial flow velocity profile develops dramatic structures at the 
QGP-MP and MP-HRG interfaces [126[|223j . This leads to large velocity gradients across these interfaces, 
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Figure 4.4: Time evolution of the average radial flow velocity (vt) = {v±) in central Au+Au collisions, 
calculated with EOS I (left panel) and SM-EOS Q (right panel). Solid (dashed) lines show results from ideal 
(viscous) fluid dynamics. The initially faster rate of increase reflects large positive shear viscous pressure 
in the transverse direction at early times. The similar rates of increase at late times indicate the gradual 
disappearance of shear viscous effects. In the right panel, the curves exhibit a plateau from 2 to 4fm/c, 
reflecting the softening of the EOS in the mixed phase. 



inducing large viscous pressures which work to reduce these gradients. In effect, shear viscosity softens 
the first-order phase transition into a smooth but rapid cross-over transition. These same viscous pressure 
gradients cause the fluid to accelerate even in the mixed phase where all thermodynamic pressure gradients 
vanish (and where the ideal fluid therefore does not generate additional flow). As a result, the lifetime of 
the mixed phase is shorter in viscous hydrodynamics, as also seen in the right panel of Figure l473l 

Fig. 14.41 plots the time evolution of radial velocity (ut) (radial flow), calculated as an average over the 
transverse plane with the Lorentz contracted energy density 7^ e as weight function, for ideal and viscous 
fluids. One flnds that radial flow builds up more quickly in the viscous fluid. The shear viscosity driven 
acceleration is strongest in the early and middles stages, but almost negligible in the later hadronic stage, 
where the ideal and viscous fluid lines are almost parallel with each other. Compared with the old radial 
flow flgure for central Cu-|-Cu collisions in Ref. [I26j (see Fig. 6 there), one flnds a slightly larger increase 
in radial flow for the smaller Cu-|-Cu system, which can be interpreted as larger viscous effects for smaller 
systems. Such system size effects are more significant for non-central collisions, which will be discussed below. 

Non-central collisions : 

We now turn to non-central collisions and take full advantage of the ability of VISII2+1 to solve the 
transverse expansion in 2 spatial dimensions. To evaluate the anisotropic evolution of the fireball, we 
will study the time evolution of flow anisotropy (Iwa;! — fFig 14. 5p and the time evolution of "spatial 
eccentricity" and "momentum anisotropy" fFig |4.6p . The spatial eccentricity [T^lIM] characterizes the spatial 
deformation of the fireball in the transverse plane; it is defined as ex — ^i^2_|_^2^ (here (...) means averaging 
over the transverse plane with the energy density e{x) as weight function). The momentum anisotropy 

/rpxx _rpyy\ 

Ep = ^ji°cx^_j?««^ (where (...) means average over the transverse plane) |72[|93j measures the anisotropy of the 
transverse momentum density due to anisotropics in the collective flow pattern, which includes only the ideal 
fluid part of the energy momentum tensor. In viscous hydrodynamics, we also define the total momentum 
anisotropy e'^ — jpp^^s^jf^, similarly defined in terms of the total energy momentum tensor T^'^ = Tq ~\-tt^'^ , 
which additionally counts anisotropic momentum contributions arising from the viscous pressure tensor. 
Since the latter quantity includes effects arising from the deviation 5f of the local distribution function from 
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Figure 4.5: Time evolution of the transverse flow anisotropy — I'fj/I) (top row) and of the anisotropy in the 
transverse source term (jiS"^^! — (bottom row). Both quantities are averaged over the transverse plane, 
with the Lorentz-contracted energy density 7j_e as weight function. The left (right) column shows results 
for EOS I (SM-EOS Q), with solid (dashed) lines representing ideal (viscous) fluid dynamical evolution. 



its thermal equilibrium form which, according to Eq. p.Sp , also affects the final hadron momentum spectrum 
and elliptic flow, it is this total momentum anisotropy that should be studied in viscous hydrodynamics if one 
wants to understand the evolution of hadron elliptic flow. In other words, in viscous hydrodynamics hadron 
elliptic flow is not simply a measure for anisotropics in the collective flow velocity pattern, but additionally 
reflects anisotropics in the local rest frame momentum distributions, arising from deviations of the local 
momentum distribution from thermal equilibrium and thus being related to the viscous pressure. 

Fig. 14.51 shows the time evolutions of the flow anisotropy — and source term anisotropy {\S^^\ 
— jiS'^^l). In central collisions these quantities vanish. In ideal hydrodynamics, the flow anisotropy is driven 
by the anisotropic gradients of the thermodynamic pressure. In viscous fluid dynamics, the source terms 
p.2b|3.2cp . whose difference is shown in the bottom row of Fig. [4751 receive additional contributions from gra- 
dients of the viscous pressure tensor which contribute their own anisotropics. Fig. l4.5l demonstrates that these 
additional anisotropics increase the driving force for anisotropic flow at very early times {t—tq <3fm/c), 
but reduce this driving force throughout the later evolution. At times t—tq >5fm/c the anisotropy of the 
effective transverse pressure even changes sign and turns negative, working to decrease the flow anisotropy. 
As a consequence of this, the buildup of the flow anisotropy stalls around r— ro~4fm/c (even earlier for 
SM-EOS Q where the flow buildup stops as soon as the flreball medium enters the mixed phase) and pro- 
ceeds to slightly decrease therafter. This happens during the crucial period where ideal fluid dynamics still 
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Figure 4.6: Time evolution for the spatial eccentricity Cx, momentum anisotropy Cp and total momentum 
anisotropy e'^ (see text for definitions), calculated for 6 = 7fm Au+Au collisions with EOS I (left column) 
and SM-EOS Q (right column). Dashed lines are for ideal hydrodynamics while the solid and dotted lines 
show results from viscous hydrodynamics. See text for discussion. 



shows strong growth of the flow anisotropy. By the time the fireball matter decouples, the average flow ve- 
locity anisotropy of viscous hydro lags about 10-15% behind the value reached during ideal fluid dynamical 
evolution. 

Figure l4!6l correlates the decrease in time of the spatial eccentricity with the buildup of the momentum 
anisotropics Cp and Cp. In viscous dynamics the spatial eccentricity is seen to decrease initially faster than 
for ideal fluids, due to the faster radial expansion shown in Fig. 14.41 The effects of early pressure gradient 
anisotropics are reflected in the initial growth rate of the flow anisotropy (Fig. 14. 5p and the momentum 
anisotropy €p (bottom panels in Fig. 14. 6p which are seen to slightly exceed that observed in the ideal fluid at 
times up to about 3fm/c after the beginning of the transverse expansion. Figure [53] also shows that in the 
viscous fluid the flow velocity anisotropy stalls about 2fm/c (for SM-EOS Q) after start and remains about 
10% below the final value reached in ideal fluid dynamics. This causes the spatial eccentricity of the viscous 
fireball to decrease more slowly at later times than that of the ideal fluid (top panels in Fig. 14. 6p which, 
at late times, features a significantly larger difference between the horizontal (x) and vertical (?/) expansion 
velocities. We note that the differences between ideal and viscous evolution were found to be significantly 
larger in the smaller Cu-|-Cu systems studied in Ref. |126j . 

It is very instructive to compare the behaviour of the flow-induced ideal-fluid contribution to the momen- 
tum anisotropy, Cp, with that of the total momentum anisotropy e^. At early times they are very different, 
with e'p being much smaller than ep (see insets in the lower panels of Fig. 14. 6p . This reflects very large neg- 
ative contributions to the anisotropy of the total energy momentum tensor from the shear viscous pressure 
whose gradients along the out-of-plane direction y strongly exceed those within the reaction plane along the 
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Figure 4.7: Mid-rapidity particle spectra for central Au+Au collisions, calculated with EOS I (left, gluons) 
and with SM-EOS Q (right, tt^, K'^ andp). The solid blue (red dashed) lines are from ideal (viscous) hydro- 
dynamics. The purple dotted lines show viscous hydrodynamic spectra that neglect the viscous correction 
5fi to the distribution function in Eq. (|3.5|) . i.e. include only the effects from the larger radial flow generated 
in viscous hydrodynamics. 



X direction. At early times this effect almost compensates for the larger in-plane gradient of the thermal 
pressure. The negative viscous pressure gradient anisotropy (which is even larger in the smaller Cu+Cu 
systems studied in jl26| ) is responsible for reducing the growth of flow anisotropics, thereby causing the 
flow-induced momentum anisotropy ep to significantly lag behind its ideal fluid value at later times. The 
negative viscous pressure anisotropics responsible for the difference between ep and slowly disappear at 
later times, since all viscous pressure components then become very small (see Fig. 14. 101 below). 

The net result of this interplay is a total momentum anisotropy in Au-|-Au collisions (i.e. a source of 
elliptic flow V2) that for a "minimally" viscous fluid with f = 4^ is 20-25% lower than for an ideal fluid, 
at all except the earliest times (where it is even smaller). The origin of this reduction changes with time: 
Initially it is dominated by strong momentum anisotropics in the local rest frame, with momenta pointing 
preferentially out-of-plane, induced by deviations from local thermal equilibrium and associated with large 
shear viscous pressure. At later times, the action of these anisotropic viscous pressure gradients integrates 
to an overall reduction in collective flow anisotropy, while the viscous pressure itself becomes small; at 
this stage, the reduction of the total momentum anisotropy is indeed mostly due to a reduced anisotropy 
in the collective flow pattern while momentum isotropy in the local fluid rest frame is approximately restored. 



4.3 Final particle spectra and elliptic flow V2 
Spectra for central collisions: 

After obtaining the freeze-out surface, we calculate the particle spectra from the generalized Cooper-Frye 
formula (j3.5p . using the AZHYDRO algorithm [246j for the integration over the freeze-out surface S. For 
calculations with EOS I which lacks the transition from massless partons to hadrons, we cannot compute any 
hadron spectra. For illustration we instead compute the spectra of hypothetical massless bosons ("gluons"). 
They can be compared with the pion spectra from SM-EOS Q which can also, to good approximation, be 
considered as massless bosons. 
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Figure 4.8: Ratio of the viscous correction SN, resulting from the non-equihbrium correction (5/, Eq. (j3.9p . 
to the distribution function at freeze-out, to the equihbrium spectrum A'^^q = dN^q/ {dydPpT) calculated from 
Eq. (|3.5p by setting Sf = 0. The gluon curves are for evolution with EOS I, the curves for 7r~, and p 
are from calculations with SM-EOS Q. The horizontal dashed line at 0.2 indicates breakdown of the viscous 
hydrodynamic approach. 

The larger radial flow generated in viscous hydrodynamics, for a fixed set of initial conditions (shown in 
Fig. 14. 4[) , leads to flatter transverse momentum spectra [2471ll32ll248j (at least at low pr where the viscous 
correction dfi to the distribution function can be neglected in (|3.5p ). This is seen in Figure |477[ by comparing 
the dotted and dashed lines. This comparison also shows that the viscous spectra lie systematically above 
the ideal ones, indicating larger final total multiplicity. This reflects the creation of entropy during the 
viscous hydrodynamic evolution. As pointed out in [247|,I132] . this requires a retuning of initial conditions 
(starting the hydrodynamic evolution later with smaller initial energy density) if one desires to fit a given set 
of experimental py-spectra. Since we here concentrate on investigating the origins and detailed mechanics 
of viscous effects in relativistic hydrodynamics, we will not explore any variations of initial conditions. All 
comparisons between ideal and viscous hydrodynamics presented here will use identical starting times tq and 
initial peak energy densities cq. The solid lines in Figure l477l show that in our calculations for pt ^ 2 GeV/c 
the effects from Sfi have an overall negative sign, leading to a reduction of the pr-spectra at large px relative 
to both the viscous spectra without Sfi and the ideal hydrodynamic spectra. This is true for all particle 
species, irrespective of the EOS used to evolve the fluid. 

However, the earlier blast-wave model estimate from Teaney |240j found that the correction is positive, 
growing quadratically with px- A later calculation by Dusling and Teaney, using causal viscous hydro- 
dynamics in the O-G formalism and a different (kinetic) freeze-out criterium to determine the decoupling 
surface, also found a (small) positive effect from Sfi on the final pion spectra, at least up to pT = 2GeV/c, 
for freeze-out around Tdcc ~ 130 MeV, turning weakly negative when their effective freeze-out temperature 
was lowered to below 100 MeV |137j . Inspecting Eqs. (|3.5l3.10p reveals that the viscous correction Sfi in 
Eqs. p.5l3.9p depends on the signs and magnitudes of the various viscous pressure tensor components along 
the freeze-out surface, weighted by the equilibrium part f^q.i of the distribution function. Its effect on the 
final pT-spectra (even its sign!) is not a priori obvious. In Appendix lA. 41 we will explore the origin of the 
discrepancy of the sign of Sf between the different groups by investigating the effects of different freeze-out 
surfaces in hydrodynamics and in the simple blast wave model. 

In Figure 14.81 we show the non-equilibrium contribution to the final hadron spectra in greater detail. 
The figure shows that the non-equilibrium effects from Sfi are largest for massless particles and, at high pT, 
decrease in magnitude with increasing particle mass. The assumption \Sf\ <^ /eq, which underlies the viscous 
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Figure 4.9: Differential elliptic flow V2{pt) for Au+Au collisions at 6 = 7fni. Left panel: Gluons from 
evolution with EOS 1. Right panel: tt~ , , and p from evolution with SM-EOS Q. Dashed lines: ideal 
hydrodynamics. Solid lines: viscous hydrodynamics. Dotted lines: viscous hydrodynamics without non- 
equilibrium distortion 6f of distribution function at freeze-out. 



hydrodynamic formalism, is seen to break down at high pT^ but to do so later for heavier hadrons than for 
lighter ones. Once the correction exceeds 0(20%) (indicated by the horizontal dashed line in Fig. 14. 8[) . the 
calculated spectra can no longer be trusted. 

In contrast to viscous hydrodynamics, ideal fluid dynamics has no intrinsic characteristic that will tell 
us when it starts to break down. Comparison of the calculated elliptic flow V2 from ideal fluid dynamics 
with the experimental data from RHIC 26j suggests that the ideal fluid picture begins to break down 
above pr — 2.5GeV/c for pious and above pr — 3GeV/c for protons. This phenomenological hierarchy of 
thresholds where viscous effects appear to become essential is qualitatively consistent with the mass hierarchy 
from viscous hydrodynamics shown in Fig. 14.81 

In the region <pt ^ 1-5 GeV/c, the interplay between jut- and pT-dependent terms in Eq. (|3.10p is 
subtle, causing sign changes of the viscous spectral correction depending on hadron mass and pT (see inset 
in Fig. 14. Sp . The fragility of the sign of the effect is also obvious from Fig. 8 in the work by Dusling and 
Teaney [137) where it is shown that in this pT region the viscous correction changes sign from positive to 
negative when freeze-out is shifted from earlier to later times (higher to lower freeze-out temperature) . Over- 
all, we agree with them that the viscous correction effects on the py-spectra are weak in this region [137j . 
We will see below that a similar statement does not hold for the elliptic flow. 

Elliptic flovif for non-central collisions 

The shear viscous effects on the suppression of the total momentum anisotropy discussed in the 
last section also reflect themselves in a suppression of the final particle elliptic flow V2^ which is shown in 
Figure Even for the "minimal" viscosity ^ = considered here one sees a very strong suppression of the 
differential elliptic flow W2(pt) from a viscous fluid (solid lines) compared to the ideal fluid (dashed lines). 
Both the viscous reduction of the collective flow anisotropy (whose effect on V2 is shown as the dotted lines) 
and the viscous contributions to the anisotropy of the local momentum distribution (embodied in the term 
5f in Eq. (|3.5|) ) play big parts in this reduction. The runs with EOS I (which is a very hard EOS) decouple 
more quickly than those with SM-EOS Q; correspondingly, the viscous pressure components are still large at 



46 




Figure 4.10: Left: Time evolution of the various components of the shear viscous pressure tensor, normahzed 
by the enthalpy and averaged in the transverse plane over the thermalized region inside the freeze-out surface 
(note that the normalizing factor e+p ^ decreases rapidly with time) , compared with different tt"'" initial 
conditions. Right: Differential elliptic flow V2{pt) for pions from & = 7fm Cu+Cu collisions with SM-EOS Q, 
compared with different tt™" initial conditions. The thick lines in both figures are for N-S initialization: 
7r'""(To) — 2770-™", while thin lines are for zero initialization: 7r'™"(ro) ~ 0. 

freeze-out, and the viscous correction 5f to the distribution function plays a bigger role. With SM-EOS Q 
the fireball doesn't freeze out until tt™" has become very small (see Fig. 14.101 below), resulting in much 
smaller corrections from 6f (difference between dashed and dotted lines in Fig. 14. 9p . On the other hand, 
due to the longer fireball lifetime the negatively anisotropic viscous pressure has more time to decelerate the 
buildup of anisotropic flow, so V2 is strongly reduced because of the much smaller flow-induced momentum 
anisotropy Cp. 

The net effect of all this is that, for Au+Au collisions and in the soft momentum region < 1-5 GeV/c, 
the viscous evolution with ^ = leads to a ^ 30% suppression of V2 for pions, in both the slope of its 
PT-dependence and its pT-hitegrated value. (Due to the flatter p^-spectra from the viscous dynamics, the 
effect in the pT-integrated V2 is not quite as large as for V2{pt) at fixed pt-) For protons, the larger radial 
flow generated in the viscous evolution also pushes the proton elliptic flow out towards larger pT (see inserts 
in the right panel of Fie: l4.9p . causing an additional suppression of ^^(Pr) at fixed pr- This effect is much 
weaker in Cu-|-Cu collisions in fl26] where, due to the small fireball size, the radial flow effects are less 
pronounced. 

4.4 Sensitivity to relaxation time and tt™" initialization 
Initialization of tt™" : 

Lacking input from a microscopic model of the pre-equilibrium stage preceding the (viscous) hydrody- 
namic one, one must supply initial conditions for the energy momentum tensor, including the viscous pressure 
tt'"". The most popular choice has been to initialize tt™" with its Navier-Stokes value, i.e. to set initially 
TT™" = 2?7(T™". Up to this point, this has also been our choice in Chapters 4.2 and 4.3. Ref. |134| advocated 
the choice tt™" = at time tq in order to minimize viscous effects and thus obtain an upper limit on 77/5 in a 
comparison with experimental data. In the present subsection we explore the sensitivity of the final spectra 
and elliptic flow to these different choices of initialization, keeping all other model parameters unchanged. 
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Figure 4.11: Time evolution of the two 
independent viscous pressure tensor 
components ir'^'^ and E = tt^^+tt^^, 
for central Au+Au collisions (solid 
lines), compared with their Navier- 
Stokes limits 2T]a^^ and 277(cr^^+cr?'f ) 
(dashed lines), for two values of the re- 
laxation time, = Srj/sT (thick red 
lines) and — l.brj/sT (thin black 
lines). All quantities are scaled by the 
thermal equilibrium enthalpy e+p and 
transversally averaged over the ther- 
malized region inside the decoupling 
surface. 



Figure H.lOl left shows the time evolution of the viscous pressure tensor, comparing the two different initial- 
izations 7r™"(ro) — 277(t"™ and 7r'""(ro) = 0. Differences are visible only at early times t—tq < 5T,r ~ 1 fm/c. 
After Ttt^ 0.2 fm/c, the initial differences have decreased by roughly a factor 1/e, and after several kinetic 
scattering times r^r the hydrodynamic evolution has apparently lost all memory how the viscous terms were 
initialized. 

Correspondingly, the final spectra and elliptic flow show very little sensitivity to the initialization of tt™", 
as seen in Fig. 14.101 right. With vanishing initial viscous pressure, viscous effects on the final flow anisotropy 
are a little weaker (dotted lines in Fig. 14.101 right), but this difference is overcompensated in the total elliptic 
flow by slightly stronger anisotropics of the local rest frame momentum distributions at freeze-out (solid 
lines in Fig. 14.101 right). For shorter kinetic relaxation times r^r, the differences resulted from different ini- 
tializations of tt"™ would be even smaller. 

Relaxation time r^r: 

While the finite relaxation time for the viscous pressure tensor in the Israel-Stewart formalism elim- 
inates problems with superluminal signal propagation in the relativistic Navier-Stokes theory, it also keeps 
the viscous pressure from ever fully approaching its Navier-Stokes limit tt™" — 2770''"". In this subsection we 
explore how far, on average, the viscous pressure evolved by VISH2-I-1 deviates from its Navier-Stokes limit, 
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Figure 4.12: Differential elliptic flow V2{pt) for tt from 6 = 7fm Au+Au collisions with SM-EOS Q, calcu- 
lated from viscous hydrodynamics with two different values for the relaxation time r^r . 

and how this changes if we reduce the relaxation time by a factor 2. (Another investigation of relaxation 
time effects and the approach towards the N-S limit can be found in Sec. 17.3.11 where we compare the "full" 
and "simplified" I-S equations.) 

In Figure 14.111 we compare, for non-central Au-|-Au collisions, the time evolution of the scaled viscous 
pressure tensor, averaged in the transverse plane over the thermalized region inside the freeeze-out surface, 
with its Navier-Stokes limit, for two values of Ttt, = 'irj / sT — r'j^^^^ / 2 and Tt^—t"^^^^ /A. For the larger 
relaxation time, the deviation from the Navier-Stokes limit is always less than 10%, and this fraction greatly 
decreases at later times. For the twice shorter relaxation time, the fractional deviation from Navier-Stokes 
decreases by somewhat more than a factor 2 and never exceeds a value of about 5%. 

Figure l4. 1 21 shows the corresponding elliptic flow V2. For the viscous fluid curves (both the flow anisotropy 
only and the full viscous hydrodynamics curves), one finds only very slight dependencies on relaxation times 

. This result is in sharp contrast to what we found in [126] for Cu-|-Cu collisions at b=7 fm, where we 
saw a much stronger dependence on the relaxation time. As already pointed out in |127| . the main reason 
behind this is that the calculations presented here and in fl26j use different versions of the I-S equations. 
The "full" I-S equation used here greatly reduces the dependence on r^r for final observables compared to the 
"simplified" I-S equation used in the early calculations in [126j . which lead to a much stronger sensitivity to 
T^r. This will be further investigated in Chap [73Tl 

4.5 Shear viscosity effects: dependence on system size and collision energies 

In Ref. [127j . we showed qualitatively that shear viscosity effects are larger for smaller systems, by 
comparing the viscous V2 suppression for Au-|-Au collisions at b=7 fm and Cu-|-Cu collisions at b=7 fm. 
Although the two systems have similar initial eccentricity Sx, the latter system is not a scaled transformation 
of the former one due to the Woods-Saxon parametrization in the Glauber initialization (see Chap. 3.2). 
To quantitatively study the dependence of shear viscosity effects on system size and collision energies, we 

2 2 

here set the initial energy density to a simple Gaussian function, e{x,y,To) = eg • exp(— -^^^ — ^). This 
construction creates an initially elliptic fireball with eccentricity £x{tq) = 0.6, independent of the parameter 
R. To study system size effects, we do comparison runs with i? = 4 fm, tq = 0.4 fm/c and i? = 6 fm 
To — 0.6 fm/c; the second system is an exact scale transformation of the first one in both space and time, by 
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Figure 4.13: Time evolution of total momentum anisotropy e'p from ideal and viscous hydrodynamics with 
different initial conditions. For viscous hydrodynamics: rj/s — 0.08 and t-^ — 

a factor 1.5. To study collision energy effects, we compare hydrodynamic evolutions starting with an initial 
energy density eg = 30 GeV/fm'^ and with eo = 60 GeV/fm'^. 

Fig. 14.131 shows the time evolution of total momentum anisotropy e'p for different initial conditions (with 
ideal EOS: EOS I). Instead of using the proper time t for the horizontal axis, we change it to the normalized 
unitless "time" (r — to)/R. For ideal fluid evolution, one observes exact scale invariance, independent of 
system size or collision energy, due to the scale invariance of the ideal hydrodynamic equations. However, 
realistic simulations for relativistic heavy ion collisions requires a freeze-out condition, a realistic EOS with 
a phase transition, and using more complicated initial conditions (e.g. the Glauber model). All of these 
introduce additional external scales (an external time scale through the freeze-out condition, an external 
energy scale through Tc in the EOS, the surface thickness in the Wood-Saxon density profile), which break 
the scale invariance of ideal hydrodynamics. This will be discussed in more detail in Chap. 5. 

Microscopically, a finite shear viscosity corresponds to a finite mean free path, which introduces an ad- 
ditional length scale, again breaking the scale invariance of the hydrodynamic system. If one examines the 
analytical solution (jA.9|l for the relativistic Navier-Stokes equations for a (O-l-l)-d longitudinally boost in- 
variant system (see Appendix lA.3.2|l . one finds a scale-breaking term proportional to This term shows 
that, for a constant rj/s, the scale breaking effect is smaller for larger systems or higher collision energies 
(corresponding to larger tq or Tq, respectively). For (2-fl)-d viscous hydrodynamics, there is no analytical 
solution. The numerical results in Fig. 14.131 clearly show scale breaking effects caused by shear viscosity in 
the viscous e'p suppression, relatively to the (universal) ideal fluid lines. For identical collision energy, one 
finds stronger suppression (i.e. larger scale breaking effects) in the smaller system (comparing red and black 
solid lines); for identical system size, one finds smaller viscous suppression for higher collision energies 
(comparing red and green solid lines). Both tendencies agree with qualitative expectations based on the 
analytical results for the (l-l-l)-d case. The system size and collision energy dependence of the viscous 
suppression translates into a corresponding dependence of the viscous V2 suppression, which will be further 
studied in Chap. 5. 
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Figure 4.14: Characteristic transverse momentum where the viscous corrections to the final pion spectrum 
become so large (> 20%) that the spectrum becomes unreliable, as a function {1/ S)dNch/dy. Stars are for 
central Au+Au collisions, open circles for semi-central Au+Au collisions at 6 = 7fm. 



4.6 Breakdown of viscous hydrodynamics at high px 

As indicated by the horizontal dashed lines in Figs. 14.81 the assumption <C |/eq| under which the 
viscous hydrodynamic framework is valid breaks down at sufficiently large transverse momenta. For a quan- 
titative assessment we assume conservatively that viscous hydrodynamic predictions are no longer reliable 
when the viscous corrections to the particle spectra exceed 20%. Fig. 14. 81 shows that the characteristic trans- 
verse momentum where this occurs depends on the particle species and increases with particle mass. To 
be specific, we here consider p^ for pions — the values for protons would be about 15% higher. 

In Fig. 14. 141 we show the breakdown momentum p^ for pions as a function of (1/ S)dNch/dy (final multi- 
plicity dNch/dy per overlap area S, the details for calculating [\/ S)dNch/dy can be found in Chap. lO)) . for 
both central and semi-central Au+Au collisions. (Larger {\/ S)dNch/dy corresponds larger collision energies. 
The initial time was fixed by taking the relation tqTo — const., following Ref. [249].) One finds that p'^ 
rises with collision energy for both central and non-central Au+Au collisions. The rise of with increasing 
e(r=0) reflects the growing fireball lifetime which leads to smaller viscous pressure components at freeze-out. 
This lifetime effect is obviously stronger for central than for non-central collisions, leading to the faster rise 
of the stars than the open circles in Fig. 14.141 For fixed {l/S)dNch/dy (which corresponds to similar initial 
entropy densities in the fireball center), we find that in central collisions (6 = 0) the validity of viscous 
hydrodynamics extends to larger values of pt than in non-central collisions (6 = 7 fm) : Viscous effects are 
more serious in non-central than in central collisions. 



4.7 Conclusions 

In this chapter, we numerically studied the shear viscous effects on the hydrodynamic evolution, final 
hadron spectra, and elliptic flow V2. Complementary to our early paper |1251ll26j , which concentrated on 
Cu+Cu collisions and used the"simplifled I-S equation", we studied Au+Au collisions, using an updated 
version of VISII2+1 that solves the "full I-S equation". 

By comparing ideal and viscous hydrodynamic runs, we explored the effects of shear viscosity for a 
"minimally" [146j viscous fluid with ^ = ^ in central and non-central Au+Au collisions, comparing the 
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evolution for two different equations of state - an ideal massless parton gas (EOS I) and an EOS with a 
semirealistic parametrization of the quark- hadron phase transition (SM-EOS Q). 

We found that shear viscosity decelerates longitudinal expansion, but accelerates the build-up of trans- 
verse flow. This slows the cooling process initially, leading to a longer lifetime for the QGP phase, but causes 
accelerated cooling at later stages by faster transverse expansion. For SM-EOS Q, we saw that viscous 
pressure gradients during the mixed phase increase the acceleration during this stage and slightly reduce its 
lifetime. They counteract large gradients of the radial velocity profile that appear in ideal fluid dynamics as 
a result of the softness of the EOS in the mixed phase, thereby de facto smoothing the assumed first-order 
phase transition of SM-EOS Q into a rapid cross-over transition. In the end the larger radial flow develop- 
ing in viscous hydrodynamics leads to flatter transverse momentum spectra of the finally emitted particles, 
while their azimuthal anisotropy and the final elliptic flow V2 in non-central heavy-ion collisions is found 
to be strongly reduced. These generic shear viscous effects are qualitatively similar for different systems at 
different collision energies, as seen by comparing the results for Au-|-Au collision shown in this chapter with 
our early results for Cu-fCu collision presented in Ref. [126j . On a more quantitative level, one finds that 
smaller systems and lower collision energies lead to larger viscous effects, as shown in Chap. 4.5. 

Our studies show that shear viscous effects are strongest during the early stage of the expansion phase 
when the longitudinal expansion rate is largest. The discussion in Chap. 4.5 suggests that TgO = -^{d ■ u) 
is the parameter that controls the strength of viscous effects. At later times the viscous corrections become 
small, although non-negligible. Small non-zero viscous pressure components along the hadronic decoupling 
surface have significant effects on the final hadron spectra that grow quadratically with transverse momentum 
and thus limit the applicability of the viscous hydrodynamic calculation to transverse momenta below 2- 
3GeV/c, depending on impact parameter, collision energy, and particle mass. Viscous effects are more 
important in peripheral than in central collisions, and larger for light than for heavy particles. Since the 
breakdown of viscous hydrodynamics is signalled by the theory itself, through the relative magnitude of the 
viscous pressure, the applicability of the theory can be checked quantitatively case by case and during each 
stage of the expansion. 

For the kinetic relaxation times r^r considered in the present work, sensitivities to the initial value of the 
viscous pressure tensor were found to be small and practically negligible. Sensitivity of the elliptic flow to 
the value of r^r was also found to be weak if the full I-S equation is used when solving viscous hydrodynamics. 

The most important finding here is the large V2 suppression caused by shear viscosity. While the viscous 
hardening of the hadron p^-spcctra can be largely absorbed by retuning the initial conditions, starting the 
transverse expansion later and with lower initial entropy density |131[|132] . this only acerbates the viscous 
effects on the elliptic flow V2 which is further reduced by the decreased fireball lifetime following the retuning. 
The reduction of the elliptic flow V2 by shear viscous effects is therefore a sensitive and robust diagnostic 
tool for shear viscosity of a fluid. This exciting finding can be used to extract the shear viscosity from 
experimental data. Similar conclusions have been obtained by several groups in the past three years, using 
independently developed viscous hydrodynamic codes, with different initial and final conditions and different 
EOS input ^125, 126. 134lll35lll37j . The extraction of the QGP viscosity from experimental data, together 
with an uncertainty analysis based on currently available viscous hydrodynamic studies, will be discussed in 
Chap. 8. 
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Chapter 5: Multiplicity Scaling of vijE 
and Shear Viscosity Effects 



5.1 Introduction 

This chapter wih continue to investigate shear viscosity effects by studying the multiplicity scaling of the 
normalized elliptic flow v^je within ideal and viscous hydrodynamics. The results compiled in this thesis 
are based on our paper [127j . The motivation for this research is provided by the well-known systematic 
comparison of Voloshin ei al. [2501125 11I252| of elliptic flow data with ideal fluid dynamical predictions which 
suggests that the elliptic flow parameter vi scaled by the initial source eccentricity e, W2/e, while strongly 
deviating from ideal hydrodynamics at low multiplicities, still scales with the final multiplicity dNch/dy per 
unit overlap area S: 

V2 1 dNcY, 

For ideal fluids the right hand side is a direct measure of the initial entropy density [253j . The scaling 
(|5.ip implies that all dependence on impact parameter, collision energy and system size can be, to good 
approximation, absorbed by simply taking into account how these control parameters change the final hadron 
multiplicity density. We will call this observation simply "multiplicity scaling of the elliptic flow", where 
"elliptic flow" is a shorthand for the eccentricity-scaled elliptic flow V2/e and "multiplicity" stands for 
{l/S)dN,y,/dy. 

Such a scaling is expected for ideal fluid dynamics whose equations of motion are scale invariant and 
where the eccentricity-scaled elliptic flow is therefore predicted [72ll254j to depend only on the squared speed 
of sound, = §f , which describes the stiffness of the equation of state (EOS) or "pushing power" of the 
hydrodynamically expanding matter. It has been known, however, for many years [223] that this ideal-fluid 
scaling is broken by the final freeze-out of the matter: if hadron freeze-out is controlled by hadronic cross 
sections (mean free paths) or simply parametrized by a critical decoupling energy density edoc or temperature 
Tdoc, this introduces an additional scale into the problem that is independent of (or at least not directly 
related to) the initial geometry of the fireball. This breaks the above argument based on scale invariance of 
the ideal fluid equations of motion. We will show here that this also leads to a breaking of the multiplicity 
scaling of V2/S not only in the most peripheral or lowest energy collisions, where freeze-out obviously cuts 
the hydrodynamic evolution short since the freeze-out density is reached before the flow anisotropy can 
fully build up |223', but even in the most central collisions at RHIC where freeze-out still terminates the 
hydrodynamic evolution before the elliptic flow can fully saturate (see also [12811229] ). 

The more interesting aspect of the experimentally observed scaling is, however, its apparent validity in 
regions where ideal fluid dynamics does not work (these encompass most of the available data [250j ). Many 
years ago, simple scaling laws for the centrality dependence of elliptic flow were derived from kinetic theory 
in the dilute gas limit, where the particles in the medium suffer at most one rescattering before decoupling 
[25511256] : these can be reinterpreted in terms of multiplicity scaling for V2/S. The dilute gas limit is expected 
to hold for very small collision systems, very large impact parameters or very low collision energies. More 
recently, a successful attempt was made to phenomenologically connect the dilute gas and hydrodynamic 
limits with a 1-parameter fit involving the Knudsen number [257] . This fit works very well for Au-f Au and 
Cu-|~Cu data from RHIC, but predicts that even in the most central Au-|-Au collisions at RHIC the ideal fluid 
dynamical limit has not yet been reached and is missed by at least 25% [257] . In this chapter, we use viscous 
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relativistic hydrodynamics to explore the multiphcity scahng of v^je in the phenomenologically relevant 
region. We conclude (not surprisingly since much of the available data is from regions where the viscous 
hadronic phase plays a large role [229]) that the multiplicity scaling data |250[|251[|252| require significant 
shear viscosity for the medium, especially during its hadronic stage, but also that viscous hydrodynamics 
predicts subtle scaling violations which seem to be qualitatively consistent with trends seen in the data (even 
if the experimental evidence for scaling violations is presently not statistically robust) and whose magnitude 
is sensitive to the specific shear viscosity 77/5. This gives hope that future, more precise data can help 
constrain the QGP shear viscosity through exactly such scaling violations. 

We should caution the reader that, similar to Ref. [257j which used a constant (time-independent) 
cross section for the Knudsen number, our viscous hydrodynamic calculations are done with a constant 
(temperature-independent) specific shear viscosity ry/s. Neither assumption is realistic, and we expect 77/5 
in particular to show strong temperature dependence near Tc (the critical temperature for the quark-hadron 
phase transition, see Chap. 1.4 for details) and emerge from the phase transition with much larger values 
than in the QGP phase. Comparisons between the results presented here and experimental data are there- 
fore, at best, indicative of qualitative trends, and improved calculations, which in particular match viscous 
hydrodynamics to a realistic hadron cascade below Tc, are required before quantitative extraction of 77/s 
from experimental data can be attempted. 

5.2 Multiplicity scaling of v^j^ in ideal and viscous hydrodynamics 

In this section we explore the multiplicity scaling (as defined in the introduction) of the eccentricity-scaled 
elliptic flow v^js, comparing ideal fluid dynamics with that of near- minimally viscous fluids with specific 
shear viscosity f = O (^) with three different EOS: EOS I, SM-EOS Q and EOS L (see FiglXHand expla- 
nation there.). Here we consider both Au+Au collisions and Cu+Cu collisions. The initial conditions, final 
conditions and other free inputs are described in Chap. 3. In the spectra calculation, we neglect resonance 
decays and only show the elliptic flow of directly emitted pions. To estimate the total charged hadron mul- 
tiplicity, we take the directly emitted positive pions, multiply by 1.5 to roughly account for multiplication 
by resonance decays at Tdcc, then multiply by another factor 2 x 1.2 = 2.4 to account for the negatives and 
roughly 20% of final charged hadrons that are not pions. A proper calculation of the resonance decay chain 
is computationally expensive and, for a systematic study like the one presented here that requires hundreds 
of runs of VISH2+1, beyond our presently available resources. 

5.2.1 EOS I: conformal fluids with e = 3p 

We begin with the simple case of a conformal fiuid with the equation of state e = 3p (EOS I), without 
phase transition. In this case the speed of sound is a constant, independent of temperature T, = For 
the ideal fluid case, naive scaling arguments based on the scale invariance of the ideal fluid equations of 
motion would thus predict a constant W2/e, independent of multiplicity density {1/ S)dNch_ldy. (The nuclear 
overlap area S is computed as S' = tt^J (x^) (y^) where (...) denotes the energy density weighted average over 
the transverse plane, and we here use the standard the eccentricity e = fa2+x2) ? which is slightly different 
from the one we used in [127j . The left panel in Fig. 15.11 clearlv contradicts this expectation. Freeze-out at 
Tdec = 130 MeV cuts the hydrodynamic evolution of the momentum anisotropy £p short before the elliptic 
flow has fully saturated. As the left panel of the flgure shows, this not only causes a strong suppression of 
W2/e at low multiplicity densities, where the time between beginning of the hydrodynamic expansion and 
freeze-out becomes very short, but it also breaks the multiplicity scaling at high multiplicity density, albeit 
more weakly. At a fixed value of {1/ S)dNch/dy, one sees larger V2/e for more central collisions initiated at 
lower collision energies (corresponding to smaller cq parameters) than for more peripheral collisions between 
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Figure 5.1: The eccentricity-scaled elliptic flow V2/£ as a function of charged multiplicity density, 
{1/ S)dNch/dy, for a conformal fluid with EOS I. Results for Cu+Cu and Au+Au collisions with two differ- 
ent initial energy densities at a variety of impact parameters, as indicated in the legend, are superimposed. 
Results from ideal fluid dynamics (a) are compared with those from viscous hydrodynamics, using the full 
Israel-Stewart equations (b). In all cases approximate, but not perfect multiplicity scaling is observed (see 
text for discussion) . In panels (b) , the ideal fluid results from the left panel are reproduced as brown symbols 
for comparison. 



the same nuclei at higher beam energies, and also for more central Cu+Cu collisions (with a rounder shape) 
than for more peripheral Au+Au collisions (with a more deformed initial shape). We find that the larger 
V2/e values in central compared to peripheral collisions can be traced directly to somewhat longer lifetimes 
of the corresponding fireballs, i.e. to the availability of more time to approach the saturation values of 
the momentum anisotropy and elliptic fiow before reaching freeze-out [128] (see Fig. IA.5l in Appendix IA.5[) . 
These freeze-out induced scaling violations in ideal fluid dynamics disappear at sufficiently high collision 
energies (i.e. large cq) where the momentum anisotropy has time to fully saturate in all collision systems 
and at all impact parameters, before freezing out (see Fig. 15.21 below). 

The right panels in Fig. 15.11 show the analogous results for a minimally viscous fluid with f = 4^ and 
kinetic relaxation time t^^ ~ Along with the suppression of U2 /e by shear viscosity, we see the appearance 
of additional scale-breaking effects. Shear viscosity breaks the multiplicity scaling of V2/e because (as shown 
in the preceding section) viscous effects are larger in smaller collision fireballs. Consequently, if we compare 
different collision systems that produce the same charged multiplicity density (1/ S)dNch/dy, we find smaller 
V2/e for Cu+Cu than for Au+Au collisions, and for peripheral Au+Au collisions at higher collision energy 
than for more central Au+Au collisions at lower collision energy. Again, the amount of elliptic flow V2/S 
generated for given (l/S)dNch/dy correlates directly with the time before freeze-out, as shown in the right 
panel of Fig. IA.5I in Appendix IA.5I 

Viscous effects also generates entropy, which increase the final charged multiplicity dN^h/dy. Comparing 
in the right panels of Fig. 15. li the brown (shaded) symbols from ideal fiuid dynamics with the colored (solid) 
symbols for viscous hydrodynamics, points corresponding to the same collision system and impact param- 
eter are seen to be shifted to the right. This enhances the scaling violations: for a given collision system, 
impact parameter and collision energy, viscosity decreases the eccentricity scaled elliptic fiow W2/e, pushing 
the corresponding point downward in the diagram, and simultaneously increases the entropy, pushing the 
corresponding point horizontally to the right. The combination of these two effects separates the curves for 
different collision systems and energies farther in viscous hydrodynamics than in ideal fiuid dynamics. 
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Figure 5.2: Same as Fig. O but for SM-EOS Q (top row) and EOS L (bottom row). For the ideal fluid 
case (a,c) an extended range of cq values up to cq = 120GeV/fm'^ was studied, in order to show that V2/S 
eventually increases again at higher collision energies [223] . 



5.2.2 Phase transition effects: EOS Q vs. EOS L 

Figure 15.21 shows the analogous results if the fluid evolves under the influence of an equation of state 
with a quark-hadron phase transition, EOS Q (top row) or EOS L (bottom row). Again approximate 
multiplicity scaling of V2/e is observed, but some scale-breaking effects are visible in both ideal and viscous 
hydrodynamics. For the equations of state with a phase transition, the scale-breaking effects are actually 
larger in the ideal than in the viscous case, i.e. in viscous hydrodynamics V2/S shows better multiplicity 
scaling than in ideal fluid dynamics! We interpret the large scale-breaking effects in the ideal fluid case as a 
complication arising from interference between the freeze-out process and the weak acceleration of matter in 
the phase transition region. This interpretation is supported by a comparison between SM-EOS Q with its 
first-order phase transition (upper left panel in Fig. l5.2p and the smooth crossover transition in EOS L (lower 
left panel): for ideal fluids, the scale-breaking effects are obviously larger for SM-EOS Q than for EOS L. 
As already observed in [125] and in Chap. 14. 2[ shear viscosity effectively smears out the phase transition 
and reduces its effect on the dynamics. In Fig. 15.21 this is clearly seen on the left side of each panel (i.e. at 
small values of ^^^^) where for the ideal fluid V2/e shows a non-monotonic peak structure [223] that is 
completely gone in the viscous case. 
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Figure 5.3: (a) The experimental observation of multiplicity scaling for W2/e, with data from Au+Au, 
Pb+Pb, and Cu+Cu collisions at RHIC, SPS and AGS 251,252 . (b) Theoretical prediction of approximate 
multiplicity scaling from viscous hydrodynamics using the full I-S equations, for three different (constant) 
specific entropy values ij/s = 0.08, 0.16, 0.24. 



It is interesting to observe that, for ideal fluids, EOS L leads to about 10% more elliptic flow under RHIC 
conditions than SM-EOS Q. The reason is that in the phase transition region EOS L is stiffcr than SM- 
EOS Q. This plays an important role at RHIC because the softness of the EOS near Tc inhibits the buildup 
of elliptic flow exactly under RHIC conditions [223j . As a corollary we note that, if RHIC elliptic flow data 
exhaust ideal fluid predictions made with SM-EOS Q [26], they will not exhaust ideal fluid predictions based 
on EOS L, thus leaving some room for shear viscous effects. 



5.2.3 A look at the experimental data 

Figure [?751 a) shows the famous experimental plot by Voloshin [251ll252j which provides empirical evidence 
for multiplicity scaling of V2/e. The lines labelled "HYDRO" are sketches for expectations from rough ideal 
fluid dynamics, based on early calculations presented in [223] for V2 in Au+Au collisions at fixed impact 
parameter = 7fm as a function of multiplicity (parametrized by eo). They should be replaced by the curves 
shown in the left panels of Fig. 15.21 

In Fig. I5.3r b) we present multiplicity scaling curves for V2/S obtained from viscous hydrodynamics with 
the full I-S equations. On a superficial level, the theoretical curves show qualitative similarity with the 
experimental data, giving correct ball-park numbers if one assumes 77/s ~ 0.24 ^ 3/47r. Interestingly, 
ignoring experimental error bars, one can see evidence for small scaling violations in the experimental data 
whose pattern agrees with the theoretical predictions from viscous hydrodynamics (see discussion at the 
end of Sec. 15.2. ip : the 62.5 A GeV Au+Au data lie shghtly above the 200 A GeV Au+Au points, and the 
200 A GeV Cu+Cu points fall slightly below the 62.5 A GeV Au+Au data. Of course, these fine features of 
the experimental data are presently not statistically robust (significant); much more precise data are needed 
to confirm or disprove the theoretical predictions, but upcoming high-statistic runs at RHIC should be able 
to deliver them. 

Closer inspection of the two panels in Fig. 15.31 shows, however, that the theoretical scaling curves have 
the wrong slope: on the left side of the plot, i.e. for small multiplicity densities, the data seem to point 
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towards larger specific shear viscosities 7 > 3 x whereas on the right side of the plot, for ^ ^^^i^ > 20 fm ^ , 
the experimental data require smaller shear viscosities, 7 < (1—2) x But this is not at all unexpected: 
Collisions represented by points in the right half of the plot correspond to high collision energies and large 
initial energy densities whose expanding fireballs spend the largest fraction of their life in the QGP phase. 
Fireballs created in collisions represented by points in the left part of the diagram have smaller initial energy 
densities and thus spend most of their time in the much more viscous hadronic phase [229j . A meaningful 
comparison between theory and experiment thus must necessarily account for the temperature dependence 
of rj/s and its dramatic increase during the quark-hadron phase transition [37 j . This would lead to scaling 
curves in Fig. IS.Sf b) with a larger slope that can better reproduce the data. What one can say already now 
is that the high-energy end of Fig. 15.31 requires very small specific shear viscosity ij/s for the QGP, of the 
same order as the minimal value postulated in |146| (unless, by using the Glauber model eccentricity instead 
of the fKLN eccentricity, the initial source eccentricity e was strongly underestimated in the experimental 
data, see discussion in Chap. 8). 



5.3 Multiplicity scaling of entropy production in viscous hydrodynamics 

In this section, we investigate the entropy production in viscous hydrodynamics for Au+Au and Cu+Cu 
systems at a variety of collision energies and impact parameters. Similar to Ref. |132j . we compute entropy 
production by exploiting the proportionality of final entropy to final charged multiplicity. We compute the 
final multiplicity dNct/dy for both ideal and viscous hydrodynamics and then equate the fractional increase 
in dNch/dy with the fractional increase in dS/dy. This ignores a small negative correction due to the viscous 
deviation of the distribution function on the freeze-out hypersurface from local equilibrium [20111212) which 
slightly reduces the entropy per finally observed particle in the viscous case. The real entropy production 
is thus slightly smaller than calculated with our prescription. However, since on the freeze-out surface the 
viscous pressure components are small [125j . this correction should be negligible. 

We checked the above procedure by also directly integrating the viscous entropy production rate d ■ 
s = n^''n^^/2r] over the space-time volume enclosed between the initial condition Cauchy surface and the 
final freeze-out surface. This method results in slightly larger entropy production, the relative difference 
amounting to about 0.7% (or about 0.07% in the absolute value of AS/Sq) for central Au-f Au collisions. 
Since the estimate from the final multiplicity gives a lower entropy production value even without accounting 
for the somewhat smaller entropy per particle in the viscous case, we conclude that entropy production due to 
numerical viscosity must be a bit smaller in the viscous fluid than in the ideal one. This is not unreasonable, 
given the observation in [125] that, compared to the ideal fluid case, the physical viscosity smoothens the 
strong velocity gradients near the quark-hadron phase transition, thereby presumably also reducing the 
effects of numerical viscosity. 

Figure [5^ shows the viscous entropy production AS, as a fraction of the initial entropy Sq, for Cu-|-Cu 
and Au-f Au collisions at various impact parameters and collision energies, as a function of multiplicity 
density. One observes approximate multiplicity scaling of the fractional entropy production, with scaling 
functions that depend on the equation of state. As for V2/e, we see small scale-breaking effects, but generally 
the produced entropy fraction shows better multiplicity scaling than elliptic flow. The scale breaking effects 
for the viscous entropy production rate go in the same direction as with elliptic flow insofar as, at the same 
value of ^ , larger collision systems and more central collisions produce fractionally more entropy than 
smaller or more peripheral collisions, due to their longer lifetimes before freeze-out (see Fig. IA.5] ). Figure [531 

We note that our viscous evolution starts earlier (at tq = OSfm/c) than that of Ref. \132f (who use tq = Ifm/c). This 
earlier start results in larger entropy production fractions. As the inset in Fig. \ 5.4}{ b) shows, most of the entropy is produced 
during the early stage of the expansion. We have confirmed that the difference between Ref. \13SH and the work here is 
quantitatively reproduced by the entropy generated during the time interval from 0.6 to l.Ofm/c, which can be calculated to 
excellent approximation analytically t211( (using Eq. (D3) in Ref. fl25\l ) by assuming boost-invariant longitudinal expansion 
without transverse flow during this period. 
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Figure 5.4: Entropy production AS, normalized by the initial entropy 5o, as a function of charged multiplicity 
density -| . Calculations with VISH2 + 1 were performed for Au+Au and Cu+Cu collisions at various 
impact parameters and collision energies, using 77/s = 0.08, r^r — Srj/sT, and three different equations of 
state (EOS I, SM-EOS Q, and EOS L). The inset in panel (b) shows the entropy production as a function 
of time, for central Au+Au collisions with parameters as indicated in the legend. 



also shows 15 — 25% differences between the entropy production rates for EOS Q (first order phase transition) 
and EOS L (rapid crossover transition). The differences are largest for the most central Au+Au and Cu+Cu 
collisions at top RHIC energies. The somewhat stiffer nature of EOS L near Tc causes the fireball to expand 
faster and with higher acceleration, leading to larger viscous effects than for EOS Q. 

An important comment relates to the negative overall slope of the scaling curves for entropy production 
shown in Figs. 15.41 Since peripheral collisions produce relatively more entropy than central collisions, and 
the produced entropy is reflected in the final charged hadron multiplicity, the collision centrality dependence 
of hadron multiplicities is altered by viscous effects. When viscous efl'ects are accounted for, the charged 
multiplicity dNch/dy will rise more slowly as a function of the number of participant nucleons A'part than 
for an ideal fluid with the same set of initial conditions. In a Glauber model parametrization of the initial 
conditions [26; this tempering effect will have to be compensated for by increasing the "hard" component in 
the initial entropy production, i.e. the component that scales with the density of binary collisions and is thus 
responsible for the non- linear increase of dNch/dy with iVpart- In the Color Glass Condensate approach [225] 
this non-linear rise is controlled by the centrality dependence of the saturation momentum scale Qs, with 
no free parameters to tune. It remains to be seen whether the success of the CGC model in describing the 
centrality dependence of dNch/dy '258] survives the inclusion of entropy (or multiplicity) producing effects 
resulting from shear viscosity during the evolution from the initial CGC to the finally observed state. 



5.4 Concluding Remarks 

The main motivation for the work presented in this chapter was provided by the experimentally observed 
multiplicity scaling of the elliptic fiow, shown in Fig. 15.3b . and its deviation at low multiplicities from ideal 
fluid dynamical predictions. We saw that many of the observed features are qualitatively consistent with 
viscous hydrodynamic calculations as presented here, and that the same calculations also predict approximate 
multiplicity scaling for viscous entropy production. Our studies revealed, however, that even for ideal fluid 
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dynamics the multiplicity scaling of the elliptic flow is not perfect, with small scaling violations introduced by 
the initial density profiles and by the freeze-out process which cuts the evolution of elliptic flow short. Even 
at RHIC energies, where the elliptic flow almost saturates before freeze-out, kinetic decoupling truncates the 
momentum anisotropy at values slightly below their asymptotic saturation value, and the deviations depend 
on the size of the colliding nuclei and the deformation of the fireball created in the collision through the time 
available for building elliptic flow before freeze-out. 

Shear viscosity strongly suppresses the build-up of momentum anisotropy and elliptic flow, especially for 
low multiplicity densities, i.e. at large impact parameters, low collision energies or for small sizes of the 
colliding nuclei. This changes the slope of the multiplicity scaling curve for vijs but preserves, to good 
approximation, its general scaling with -j^^^*^- Violations of multiplicity scaling for W2/e are somewhat 
larger for the viscous expansion than for the ideal fluid (especially with EOS I), but remain small enough to 
be consistent, within statistical errors, with the experimental observation of approximate scaling. The slope 
of the approximate scaling curve and the spread around this curve caused by scaling violation increase with 
the value of the specific shear viscosity 77/5 and can thus be used to constrain it. 

Specifically, the observed scaling violations have the following features: At fixed multiplicity density 
viscous hydrodynamics predicts slightly larger elliptic flow V2IS for larger collision systems or more 
central collisions than for smaller nuclei colliding at similar energy or more peripheral collisions between 
similar-size nuclei colliding at higher energy. Larger vije values are associated with longer lifetimes of the 
corresponding fireballs before freeze-out and thus also with larger relative entropy production. This correlates 
the scaling violations for vije observed in Figs. 15.21 and 15.31 with those for the relative entropy production 
AiS/iSo seen in Fig. 15.41 The pattern of the predicted scaling violations shows qualitative agreement with 
experiment, although higher quality data are required to render this agreement statistically robust and 
quantitative. 

For a fixed (i.e. temperature independent) ratio 77/s, the slope of the multiplicity scaling curve for 
Vije does not agree with experiment - the curves predicted by viscous hydrodynamics are too flat. The 
slope can be increased by allowing r\l s to increase at lower temperatures: For small multiplicity densities 
(very peripheral collisions or low collision energies), the data seem to require ^ > 3 x whereas at large 
multiphcity densities they appear to constrain the specific shear viscosity to values of ^ < (1—2) x While 
this is qualitatively consistent with the idea that in high-multiplicity events the dynamics is dominated by 
the QGP phase (whose viscosity would thus have to be small, of order l/47r) whereas low- multiplicity events 
are predominantly controlled by hadron gas dynamics (which is highly viscous [229j ). much additional work 
is needed to turn this observation into quantitative constraints for the function ^(r). 
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Chapter 6: Generic Viscous Effects 
— Bulk Viscosity 



6.1 Introduction 

In Chap. 4, we studied the generic shear viscosity effects and showed that eUiptic flow V2 is very sensitive 
to the shear viscosity to entropy ratio rj/s, and can thus be used to extract the QGP shear viscosity from 
experimental data |125|ll26|.ll27|.ll35|ll36|ll37j . In this chapter, we will investigate the bulk viscosity effects 
and discuss the uncertainties introduced by bulk viscosity when trying to constrain the QGP shear viscosity 
from experimental data |1291ll30j . 

We will first compare the generic shear and bulk viscosity effects on the hydrodynamic evolution and final 
particle spectra and elliptic flow, by comparison runs from a) ideal hydrodynamics, b) viscous hydrodynamics 
with only shear viscosity and c) viscous hydrodynamics with only bulk viscosity. In these comparison runs, 
identical initial conditions and final conditions are used, which are the same as in Chap. 4 and explained 
in Chap. 3. For viscous hydrodynamics, the transport coefficient are set as follows: b): rj/s = l/iir, 

= 3ry/sT = 3/47rr and C/s = (C = 0); c): = 0, C = 1 for C/s as shown in Fig. [3l3l rn = 3/47rr. 
We here use N-S initializations for both shear stress pressure tt™" and bulk pressure 11 here. In contrast 
to Chap. 4, which used EOS I (ideal EOS) and SM-EOS Q (a quasi-first order phase transition EOS), we 
employ here the lattice-based EOS L, for consistency with the bulk viscosity, which is constructed using the 
squared speed of sound from EOS L. 

Chap. 6.3 concentrates on bulk viscosity effects and studies their sensitivity to the relaxation time and 
initialization for the bulk pressure. Chap. 6.4 investigates the upper limit of C/s, consistent with the validity 
of 2nd order viscous hydrodynamics. In Chap. 6.5 we will study how much uncertainty bulk viscous effects 
might contribute to the extracted value of shear viscosity from data. 

6.2 Shear viscosity vs. bulk viscosity effects 

6.2.1 Hydrodynamic evolution 

The left panel of Fig l6.1l shows the time evolution of local temperatures from ideal and viscous hydrody- 
namics. When compared with the ideal fluid, the shear viscosity slows down the cooling process during the 
early stage, but speeds it up during the middle and late stages. As discussed in Ref. |125II126| and ChaD. [4?2l 
this is caused by the competition between decelerated longitudinal expansion and accelerated transverse 
expansion. During the early stage and before the transverse flow fully builds up, the system experiences 
quasi one-dimensional expansion, and the reduced work along the longitudinal direction (caused by shear 
viscosity) slows down the cooling process. In its middle and late stages, the system is in full 3-dimensional 
expansion, and the faster transverse expansion caused by the additional work done by shear viscous forces in 
the transverse directions leads to faster cooling. In contrast, bulk viscosity decelerates both the longitudinal 
and transverse expansion during the entire time evolution, since the negative bulk pressure effectively softens 
the EoS, especially near the phase transition. As a result the cooling process is slowed down during the whole 
evolution, with most of the effects concentrated around the phase transition where C/s peaks. 

The shear-viscosity-accelerated and bulk-viscosity-decelerated transverse expansion also manifests itself 
in the evolution of the radial flow, shown in the right panel of Fig l6.ll Compared with the ideal fluid, 
minimal shear viscosity leads to ~ 5% more radial flow (vr) at the end of the evolution. The increase of (vr) 
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Figure 6.1: Time evolution of local temperature (left) and average radial flow (right), from ideal hydrody- 
namics, viscous hydrodynamics with only minimal shear or bulk viscosity, calculated with identical initial 
and final conditions. 



mostly happens early, when the shear viscosity effects are largest. In contrast, bulk viscosity counteracts the 
build-up of radial flow, and minimal bulk viscosity reduces the final radial flow by 5% relative to the ideal 
fluid case. Most of the bulk viscous effects occur during the middle stage from 5-10 fm/c, which is when a 
large portion of the matter passes through the phase transition region where C/s is large. 

We now turn to non-central collisions. To describe the fireball deformation in configuration and momen- 
tum space, we use the spatial eccentricity ex and the momentum anisotropics Sp, e'p, following Ref. |223pi26] 
and our discussions in Chap. 4. Figure I6.2r a) shows the time evolution of the spatial eccentricity ex for 
non-central Au-f Au collisions at b=7 fm. Compared with the ideal fluid, bulk viscosity slows down the 
decrease of the spatial eccentricity Sx since the negative bulk pressure decelerates both the longitudinal and 
transverse expansion. In contrast, the spatial eccentricity of the shear viscous fluid without bulk viscosity 
drops initially faster than that of the ideal fluid due to the faster radial expansion shown in Fig. 16.11 (right), 
but decreases more slowly than the ideal fluid at later times, due to the smaller flow anisotropy created in 
the viscous dynamics (see Fig. 14.51 in Chap. 4). This is also reflected in the late-time momentum anisotropy 
shown in Fig. I6.2r b). as one sees by comparing the dashed blue ideal fluid line with the dash-dotted red 
shear viscous fluid line. 

It is worth discussing the time evolution of the momentum anisotropics ep and e'p in more detail. In vis- 
cous hydrodynamics, the difference between the total momentum anisotropy e'p (solid lines) and flow induced 
momentum anisotropy ep (dashed dotted lines) reflects the viscous pressure anisotropy. At freeze-out, this 
difference generates additional viscous effects on the V2 suppression, arising from the non-equilibrium part 
of the distribution function 6f along the freeze-out surface, which adds to or subtracts from that caused by 
the viscous suppression of the flow anisotropy. Shear viscous effects are mostly concentrated at earlier times, 
where they greatly reduce e'p in comparison with ep. The negative momentum anisotropy from the shear 
stress pressure slows down the growth of the flow anisotropy and causes the total momentum anisotropy ep 
(solid red line) to fall behind that of the ideal fluid (dashed blue line) at late times. Similarly bulk viscosity 
also suppresses the development flow anisotropy, since the gradients of the effective pressure p -I- 11 are reduced 
by the negative bulk pressure contribution. As a result, the momentum anisotropy ep from the bulk viscous 
fluid (dashed dotted brown line) stays always below the ideal hydro one (dashed blue line) . Most of the bulk 
viscous effects happen during the middle stage where e'p is seen to signiflcantly deviate from ep (comparing 
the solid and dot-dashed brown lines). In contrast to the shear viscous fluid, for the bulk viscous fluid e' 
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is larger than Cp; This is mostly due to the negative bulk pressure contribution in the denominator of the 
expression for e^. As a result, Cp from the bulk viscous fluid is even slightly larger than in the ideal fluid case. 
But this will not lead to an increase of hadronic V2 (even if one freezes the system out already at the phase 
transition), since most of the enhancement happens in the inner core of the QGP phase at temperatures 
T > 190 MeV, as the symbols along the bulk viscous fluid lines indicate. At later times, the total e'p from 
the bulk viscous fluid is always below the ideal fluid one, indicating a bulk viscous suppression of hadronic 
elliptic flow V2 after freeze-out. For Tdec = 130 MeV, all of the V2 suppression comes from flow effects, since 
in our construction the bulk viscosity ( is effectively zero below T < 150 MeV . If the system freezes-out 
near the phase transition around 165 — 175 MeV, the difference between e'p and Cp along the bulk viscous 
fluid lines indicates that the viscous corrections Sf along the freeze-out surface will increase V2', this was 
also pointed out by Monnai and Hirano |241j . who calculate Sf in terms of 11 and tt'"" along the freeze-out 
surface using Grad's 14-moment method. However, their calculation did not include the flow effects which 
(according to Fig l6.2p for Tdec = 165 — 175 MeV are least 2-3 times larger than those from Sf, leading to a 
final overall suppression of V2 ■ 



6.2.2 Spectra and elliptic flow 

The left panel in Fig. 16.31 left shows the pion pt spectra from central Au-|-Au collisions, calculated from 
ideal hydrodynamics, viscous hydrodynamics with only minimal shear viscosity, or viscous hydrodynamics 
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Figure 6.3: pr spectra and V2 for directly emitted pions (ie. without resonance decay contribution) 



with only minimal bulk viscosity, using identical initial and final conditions. We find that shear viscosity 
results in flatter spectra due to larger radial flow, as shown in Fig. 16.11 (right). In contrast, bulk viscosity 
leads to steeper spectra since it decreases radial flow build-up during the fireball expansion. At low pT, both 
shear and bulk viscosity lead to a larger yield, when compared with the ones from ideal hydrodynamics. 
This corresponds to an increase of the final multiplicity caused by viscous entropy production. Even though 
the shear and bulk viscosity act against each other in the slope of the final spectra, their effects add in the 
viscous entropy production, and an adjustment of initial conditions from ideal to viscous hydrodynamics 
can thus not be avoided if one intents to describe a given set of experimental data. In Ref. [133] . it was 
shown that by increasing the thermalization time tq and adusting the initial energy or entropy density, one 
can reduce the initial total entropy and reduce the buildup of the radial flow, thereby compensating for the 
viscous entropy production and allowing for a successful fit of the experimental measured final multiplicities 
and spectral slopes. 

In the right panel of Fig. 16. 3i we show the px dependent V2 of pions for non-central Au-I-Au collisions at 
b=7 fm. The insert enlarges the low-px region px < 0.5 GeV. We find that both shear and bulk viscosity 
suppress V2 at low px- At higher px (~ 2 GeV), shear viscosity further suppresses V2 due to a negative 
contribution Sf along the freeze-out surface. In contrast, bulk viscosity increases V2 above 1 GeV, due to 
its steepening effects on the px spectra there. At px — 0.5 GeV (approximately the mean average trans- 
verse momentum for pions) we find that minimal bulk viscosity suppresses V2 by ~ 5% and minimal shear 
viscosity suppresses V2 by ~ 20%. Using V2 data to determine rj/s, the additional 5% V2 suppression 
from minimal bulk viscosity can thus leads to ~ 25% uncertainties in the extracted value of rj/s. Therefore, 
we should not neglect bulk viscosity effects when extracting the QGP shear viscosity from experimental data. 



6.3 Bulk viscosity effects: sensitivity to relaxation time and bulk pressure ini- 
tialization 

We now concentrate on bulk viscosity and investigate the effects from different relaxation times rn and 
initializations of the bulk pressure 11. For rn we here consider three choices: the constant values rn — 0.5 
and 5fm/c, and the temperature dependent function Tn(r) ~ max[f • |(r), 0.1 fm/c] with f — 120fm/c. 
The last choice implements phenomenologically the concept of critical slowing down near a phase transition; 
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Figure 6.4: Left: Differential pion elliptic flow V2{pt) from ideal and viscous hydrodynamics, including only 
bulk viscosity. Right: Time evolution of the bulk pressure (11) averaged over the transverse plane (weighted 
by the energy density) from viscous hydrodynamics. Different curves correspond to different initializations 
and relaxation times, as indicated (see text for discussion). 



it yields rn ~ 0.6fm/c at T = 350 MeV and Tn ~ 5fm/c at T^. To study memory effects, we explore two 
different initializations for the bulk viscous pressure: (a) Navier-Stokes (N-S) initialization, H(to) = ~C,d-u, 
and (b) zero initialization, n(To) = 0. 

The left panel of Fig. 16.41 shows the differential elliptic flow V2{pt) of directly emitted pions (without 
resonance decays) for non-central Au+Au collisions at fe = 7fm, calculated from ideal hydrodynamics and 
minimally bulk viscous hydrodynamics with identical initial and final conditions. The different lines from 
viscous hydrodynamics correspond to different relaxation times rn and different initializations n(ro). One 
sees that these different inputs can lead to large uncertainties for the bulk viscous V2 suppression. For 
minimal bulk viscosity, the V2 suppression at pt = 0.5 GeV ranges from ^ 2% to ^ 10% relative to ideal 
hydrodynamics (blue dashed line in the left panel). 

For the shorter relaxation time, rn = 0.5 fmc/c, the bulk viscous V2 suppression is insensitive to the 
initialization of H, and both N-S and zero initializations show ~ 8% V2 suppression relative to ideal fluids. 
The reason behind this becomes apparent in the right panel showing the time evolution of the average 
bulk pressure (11). For short relaxation times, (11) quickly loses all memory of its initial value, relaxing 
in both cases to the same trajectory after about 1 — 2fm/c (i.e. after a few times rn). This is similar 
to what we found for shear viscosity where the microscopic relaxation times are better known and short 
{TTr{Tc) ~ 0.2 — 0.5fm/c) and where the shear pressure tensor tt™" therefore also loses memory of its 
initialization after about Ifm/c |126| . 

This changes if one accounts for the critical slowing down of the evolution of 11 near Tc. If one simply 
multiplies the constant relaxation time by a factor 10, setting rn = 5fm/c, one obtains the dotted and 
solid magenta lines in Fig. 16.41 Now the bulk viscous V2 suppression relative to the ideal fluid becomes very 
sensitive to the initialization of the bulk viscous pressure: For zero initialization n(ro) = 0, the viscous 
V2 suppression is very small (only ~ 2% at pr = 0.5GeV/c). The right panel shows that in this case the 
magnitude of the (transversally averaged) bulk pressure evolves very slowly and always stays small, leading 
to almost ideal fluid evolution. On the other hand, if 11 is initialized with its Navier-Stokes value, which 
initially is large due to the strong longitudinal expansion, it decays initially more slowly than for the shorter 
relaxation time. Its braking effect on the flow evolution is therefore bigger, resulting in much stronger 
suppression of V2 than for zero initialization, at ~ 10% slightly exceeding even the viscous V2 suppression 
seen for the tenfold shorter relaxation time. 
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Figure 6.5: Upper limits of ^/s, and entropy production as a function of relaxation time, for zero initialization 
n(To) = and N-S initialization n(To) — —(d ■ u, with varieties of conditions. 



The "critical slowing down" scenario with temperature-dependent TuiT) (black lines) interpolates be- 
tween the short and long relaxation times. As for the fixed larger value Tn = 5fm/c, V2 depends sensitively 
on the initialization of 11, but for N-S initialization the viscous V2 suppression is somewhat smaller than for 
both short and long fixed relaxation times. The reasons for this are subtle since now, at early times, the 
bulk viscous pressure 11 evolves on very different time scales in the dense core and dilute edge regions of the 
fireball. As a result, for N-S initialization the average value (11) is smaller in magnitude than for both short 
and long fixed rn, throughout the fireball evolution (right panel, black lines). 

6.4 Validity of I-S viscous hydrodynamics and upper limits for (/s 

Viscous hydrodynamic codes can not run with arbitrarily large bulk viscosity. Physically, a large enough 
bulk viscosity can lead to a negative value for the total effective pressure p -I- 11 (thermal pressure + bulk 
viscous pressure) if the expansion rate d ■ u is large. When this happens, the fluid becomes mechanically 
unstable [259j and will tend to break up. In numerical simulations, this corresponds to exponential amplifi- 
cation of local numerical errors, which will eventually stop the code from running. On the other hand, 2nd 
order viscous hydrodynamics applies only for a near-equilibrium system: 6f <^ fo- This requires that the 
dissipative flows are (much) smaller than the equilibrium energy momentum tensor. 

In this section, we will investigate upper limits of C/s allowed by 2nd order viscous hydrodynamics 
under dynamical conditions encountered in heavy-ion collisions. If the input ^/s is small, the condition 
p + n(a;, j/,r) > is stratified everywhere in the transverse {x,y) plane during the entire hydrodynamic 
evolution. As one increases (/s (which corresponds to increasing the coefficient C in our parametrization 
of {(/s){T), see Chap. 3.5), the condition p + 11 > will be violated near the phase transition in certain 
regions of {x,y) and t. In our calculations we define the upper limit of C/s as the critical value of C that 
can generate negative effective pressure, p -|- 11 < 0, at some position (x, y) and some time t, i.e. anywhere 
inside the freeze-out surface. 

Fig. 16.51 shows the upper fimit {C/s)max{Tc) as a function of relaxation time rn. The behavior of 
iC/ ^)max{Tc) as a function of rn is seen to depend strongly on the initialization of the bulk pressure. For 
N-S initialization, TI{tq) = ~Cd ■ u, {C,/ s)max{T(.) is insensitive to the relaxation time rn. In this case, the 
magnitude of the average bulk pressure 11 decreases more or less monotonically with time, as shown in the 
right panel of Fig. [6^ Violations of the positivity condition p + 11 > always happen at the starting time tq 
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and at transverse positions where the matter is near the phase transition. This leads to a {(/s)max{Tc) that 
is controlled by initial conditions and independent of the relaxation time. If one includes shear viscosity, 
iC/ s)maxiTc) depends only on the initialization time tq, but not on the value of r//s. The dependence on tq 
arises from the strong dependence of the initial bulk pressure H = —(^d ■ u on tq, through the expansion rate 
d ■ u = I/tq. Hence the maximal allowed {C,/ s)max{Tc) increases when the hydrodynamic evolution starts 
later. This tendency is illustrated by the solid red, dashed magenta and dotted orange lines in Fig. 16.51 as 
one increase tq from 0.6 fm/c to 1 fm/c and 2 fm/c, the maximal {C,l s)max{Tc) lines increases from 0.05 to 
0.09 and 0.18. 

For zero initialization n(ro) = 0, one finds a qualitatively similar dependence of (C/s)maa:(T'c) on tq: The 
curves {(,/ s)max{Tc) move up monotonically as one increases the starting time To from 0.6 fm/c to 1.0 fm/c 
(solid green and dashed dark green). In contrast to the N-S initialization, the {(,/ s)max{Tc) curves now show 
a strong dependence on relaxation time r^, rising monotonically with Ttt. The reason is that it takes some 
time for the bulk pressure 11 to develop large enough magnitudes to violate the positivity condition p+II > 0; 
again this happens typically in regions where the matter is close to the phase transition. Smaller relaxation 
times allow 11 to develop faster, resulting in a monotonic increase of {C,/ s)max{Tc) with rn- For Tn < 1 fm/c, 
we find "universal" {C/ s)max{Tc) — Tn curves that do not depend on the shear viscosity rj/s (solid black, 
green and blue curves), but move upwards as we increase the starting time tq. This is because the violation 
of the positivity condition p + 11 > generally happens at early times t < 3 fm/c when the flow profiles are 
not yet significantly affected by shear viscous effects. For the two viscous fluid lines with rj/s = 0.08 and 
0.16 (solid blue and solid green lines) one see that they continue to overlap each other for T^r > 1 fm/c, where 
they break away from the rj/s — line. For rj/s — 0, the phase transition generates larger velocity gradients 
near the phase transition, which cause instability at lower value of C/s. Shear viscosity smoothes out these 
large gradients, as discussed in Chap. 4, allowing the fluid to evolve stably up to larger values of (/s. Bulk 
viscosity C alone has no smoothing influence on sharp structures generated by a phase transition; For zero 
initialization for 11, shear viscosity is thus essential to stabilize the evolution of the viscous fluid against me- 
chanical instabilities caused by strongly negative bulk viscous pressure, especially for large relaxation time Tn. 

6.5 Extracting rj/s from experimental data: uncertainties introduced by bulk 
viscosity 

Similar to shear viscosity, bulk viscosity suppresses elliptic flow V2, as shown in Fig. 16.31 (right). The 
bulk viscosity effects, once included, will thus reduce the rj/s value needed to account for the suppression of 
a given measured amount of V2 /s below the ideal fluid expectation. Current theoretical uncertainties in the 
relaxation time Tn and the initial value of the bulk viscous pressure 11, as well as the unknown value of the 
bulk viscosity itself, thus introduce a signiflcant uncertainty in the extraction of rj/s from elliptic flow data, 
which can not be reduced without additional theoretical efforts. This is illustrated in Fig. 16. 6i where we 
investigate the additional viscous V2 suppression from minimal and maximal bulk viscosity. (Here, "maximal 
bulk viscosity" denotes the largest value compatible with p + H > everywhere, as discussed in Chap. 6.4). 
For N-S initialization (Fig. 16. 6| left), the viscous fluid lines with both shear and bulk viscosities (77/s = l/47r, 
C=1.0 or 1.3) lead to 20-25% more viscous V2 suppression at pt = 0.5 GeV than a fluid with shear viscosity 
only. Neglecting the possible need for re-tuning the initial conditions and other hydrodynamic parameters, 
this additional V2 suppression will translate into a reduction of the extracted rj/s by 0(20%), when fltting 
the same elliptic flow data. For zero initialization (right panel), the additional bulk viscous V2 suppression 
ranges from 2% to 70%. Similarly, this will reduce the the extracted ij/s by 0(2-70%). This uncertainty 
bound covers only the C/s range accessible with our viscous hydrodynamics simulations. Larger C/s value 
may be realized in nature but would completely destroy our ability to extract 77/s from a comparison of 
V2 data with viscous hydrodynamics because the later breaks down. To remove or signiflcantly reduce this 
uncertainty, it is mandatory to place strict theoretical limits on (^/s)(r) and Tn(T) as well as on the initial 
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Figure 6.6: V2(j>t) from ideal hydrodynamics, viscous hydrodynamics with minimal shear viscosity, and 
viscous hydrodynamics with minimal shear viscosity and minimal or maximal bulk viscosity. Left: N-S 
initiahzation, 7r'""(To) = 2ria""\ n(To) = — ' Light: zero initiahzation, 7r™"(To) = 0, n(ro) — 0. 



value for the bulk pressure H. 
6.6 Summary 

In this chapter, we first compared shear and bulk viscous effects by studying, for central and non-central 
Au+Au collisions, simulations based on (a) ideal hydrodynamics, (b) viscous hydrodynamics with only 
minimal shear viscosity {rj/s ~ 1/A-k) and (c) viscous hydrodynamics with only minimal bulk viscosity 
{C,/s)min- In all comparisons, we used identical initial and final conditions and the same EOS (EOS L). 
Similarities and differences between shear and bulk viscous effects are summarized as follows: Shear viscosity 
decelerates longitudinal expansion, but accelerates transverse expansion. When compared with the ideal fluid 
case, this leads to an initial slowing down of the cooling process, but to faster cooling during later stages, and 
more radial flow generation. Bulk viscosity decelerates both longitudinal and transverse expansion, resulting 
in slower cooling throughout and less creation of radial flow. Both shear viscosity and bulk viscosity generate 
entropy during the hydrodynamic evolution, and thus lead to larger final multiplicities. The px spectra, 
however, are flatter for a shear viscous fluid, but steeper for a bulk viscous fluid. For non-central Au-|-Au 
collisions, we found that both shear and bulk viscosity inhibit the development of momentum anisotropics. 
Most of the shear viscous effects happen during the early stage, while bulk viscous effects are stronger during 
the middle stage when most of the matter passes through the phase transition. As a result, both shear and 
bulk viscosity suppress elliptic flow V2 at low pT- At higher py, ^2 is further suppressed in a shear viscous 
fluid due to negative contributions from the non-equilibrium part 5f of the freeze-out distribution function. 
In contrast, V2 is slightly increased for the bulk viscous fluid above > 1 GeV, due to the steeper pT 
spectra. 

Effects arising from different relaxation times rn and different initializations of the bulk pressure 11 on the 
viscous V2 suppression were studied in Chap. 6.3, in order to assess phenomenological uncertainties arising 
from limited theoretical knowledge in the bulk viscous sector. The well-known phenomenon of critical slowing 
down during a phase transition suggests that rn may become large near Tc- If one uses larger relaxation 
times of the order of = 5 fm/c (near Tc or at all temperatures), the bulk viscous V2 suppression becomes 
sensitive to the initialization of 11, with effects ranging from ~ 2% (zero initialization for 11) to ~ 10% (N-S 
initialization for 11) for minimal bulk viscosity. This strongly suggests a need for more theoretical research 
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into the relaxation time rn and into theories that can provide reliable initializations for the bulk viscous 
pressure. 

The bulk viscosity C itself also has theoretical uncertainties: the peak value of (/s near the phase 
transition from lattice QCD estimates is around 10 times larger than the minimal AdS/CFT values, and 
even less is known for the bulk viscosity of the hadronic phase. Considering those theoretical uncertainties, 
we investigated the maximal peak value of (/s the near phase transition, (C/s)maa;(Tc), that is allowed by 
the the numerical applicability of 2nd order viscous hydrodynamics with dynamical conditions encountered 
in heavy ion collisions. For N-S initialization, {C/s)max{Tc) can be as low as 0.05 / 0.09 / 0.18, depending 
on the starting time tq for hydrodynamic stage, tq = 0.6 / 1.0 / 2.0 fm/c, and this limit is insensitive to the 
relaxation time rn. For zero initialization, {C,/ s)max{Tc) increases monotonically with the relaxation time 
Tn, and can roach values up to 0.7 at rn = 10 fm/c, for ro = 0.6 — 1.0 fm/c. 

The additional bulk viscous effects on the suppression of elliptic flow V2 (comparing to ideal hydrodynam- 
ics) were studied in Chap. 6.5. We concluded that bulk viscous effects can not be ignored when extracting 
the QGP shear viscosity to entropy ratio ry/s from experimental elliptic flow data. Pinning down the bulk 
viscous ofFects is important, especially if 77/s turns out to bo small of order of the KSS bound. It urgently 
requires to much more tightly constrain the allowed ranges for ^/s, r^r and n(ro), than presently possible. 
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Chapter 7: Recent Developments in Viscous Hydrodynamics 



7.1 Introduction 

During the past couple of years, several groups have independently developed (2+l)-dimcn- sional viscous 
hydrodynamic codes for relativistic heavy ion collisions and published their first results on shear viscosity 
effects |1251ll26lll271ll34[|135[|137] . In this chapter, I will put my own research results in context with related 
developments elsewhere, discuss similarities and differences as well as the origin of the latter. 

All groups that did viscous hydrodynamics simulations found that the elliptic flow coefficient V2 is very 
sensitive to shear viscosity: given the large expansion rates of heavy-ion collision fireballs, even minimal 
viscosity saturating the KSS bound rj/s > I/Att |146j for the viscosity to entropy density ratio can lead to a 
strong (and thus easily measurable) suppression of V2- Assuming the availability of a well-established ideal 
fluid dynamical baseline for V2 as a function of collision energy, centrality and system size, measurements of 
this suppression could thus be used to constrain the QGP shear viscosity from experimental data. 

However, the first results for viscous V2 suppression published by the different groups seemed to show 
large quantitative discrepancies, ranging from 20% to 70% even for 'minimal viscosity' rj/s = 1/An. Detailed 
inspection of the early papers of each group revealed that different initial conditions, final conditions, EOS 
and different versions of the 2nd order viscous equations were used in these calculations, and that collision 
systems of different size (Cu-|-Cu vs. Au+Au collisions) were studied. This motivated us to systematically 
investigate the physical effects of each of these differences and to collaborate with other groups to carefully 
verify our numerical codes to eliminate the possibility of numerical error. 

In this chapter, we will summarize recent progress in viscous hydrodynamics. We will show how the ap- 
parent discrepancies on viscous V2 suppression are resolved, when one properly accounts for physical effects 
arising from system size, equation of states (EOS), and different versions of the Israel-Stewart equations 
(Chap. 7.2). I will also report on the most recent code verification results, obtained within the TECHQM 
collaboration [260j (Chap. 7.4). In Chap. 7.3, we will numerically compare the different 2nd order I-S equa- 
tions obtained from the different theoretical approaches described in Chap. 2, which includes a comparison 
between the "simplified" and "full" I-S equations as well as some effects from other higher order terms. In 
Chap. 7.4, we will also numerically compare the I-S equation with the 0-G equation. Throughout this 
Chapter, we concentrated on shear viscous effects and set C = 0. 

7.2 System size, EOS and different I-S equations 

In this section, we will briefly discuss the different manifestations of shear viscosity when one varies system 
size and EOS and uses different versions of the I-S equations |127| . As mentioned in the Introduction, this 
analysis resolves the initially puzzling differences between the results published by different groups. 

Fig. 17.11 shows the differential elliptic flow V2{pt) for directly emitted pions from ideal and viscous 
hydrodynamics. Panels (a) and (b) compare two systems of different size (Cu-|-Cu and Au-|-Au, both at 
h = 7fm), using identical equations of state (SM-EOS Q), I-S equations (simphfied, eq. (I2.30|) in Chap. 2) and 
other free inputs. Although both systems have similar initial eccentricities, the smaller Cu-|-Cu system shows 
a much larger viscous V2 suppression (by almost 70% below the ideal fiuid value at = 2 GeV/c [12511126] ) 
than observed in the larger Au-|-Au system where the suppression is almost a factor two smaller. Panels 
(b) and (c) compare the same Au-|-Au system at = 7fm for two different EOS and different I-S equations. 
Changing the EOS from SM-EOS Q to EOS L reduces the viscous suppression of elliptic fiow by another 
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Figure 7.1: Differential elliptic flow V2{pt) for directly emitted pions (i.e. without resonance decay contri- 
butions), comparing results for different collisions systems and equations of state, (a) Cu+Cu at 6 = 7fm 
with SM-EOS Q. (b) Au+Au at 6 = 7fm with SM-EOS Q. (c) Au+Au at 6 = 7fm with EOS L. Dashed 
(solid) lines correspond to ideal (viscous) fluid dynamics, with parameters as indicated. 



quarter (from ~ 40% to ^ 30% below the ideal fluid value at pr — 2GeV/c). Replacing the simplified I-S 
equations used in [1251ll26j by the full I-S equations (eq. (|2.3ip in Chap. 2) used in [134] further reduces the 
V2 suppression from 30% to 25% below the ideal hydrodynamics value at pr — 2GeV/c. This final result 
is consistent with |134j . Although for EOS L the additional terms in the full I-S equations only result in 
a 5-10% difference in V2 suppression, its effect is much larger for more rapidly expanding systems, such as 
Cu-|-Cu collisions driven by a stiff conformal EOS e — 3p (see Chap. 7.4.1 and Ref. [127J for details). 

We conclude that the biggest contribution to the large difference between the results reported in Refs. |125j 
and 1 134] arises from the different collision systems studied, with much larger viscous effects seen in the 
smaller Cu+Cu system than in Au+Au collisions. The next most important sensitivity is to the EOS; for 
the most realistic EOS studied here, EOS L, the differences between using the full or simplified I-S equations 
with — 3ri/sT are only about 10% on a relative scale, or about 3% on the absolute scale set by the elliptic 
flow from ideal fluid dynamics. For smaller r^r, this last difference would shrink even further. 



7.3 Detailed comparison of different I-S equations 

As discussed in Chap. 2, the I-S formalism can be obtained from different approaches: from the 2nd 
law of thermodynamics (Chap. 2.3), from kinetic theory (Chap. 2.4) and from the conformal symmetry 
constraint (Chap. 2.5). The common terms, obtained by all of these approaches, form the so-called "simpli- 
fied" I-S equation shown in eq. (I2.30|) . The "full" I-S equation used here and in Ref. [127^ is defined as the 
simplified I-S equation plus a term that helps to maintain the conformal symmetry for a conformal fluid, 
shown in eq. ()2.31|) . With the identity (I2.32|) . which holds in the conformal limit, the "full" I-S equation can 
be written in different forms, and these variations are what one originally obtains from different 2nd order 
theory approaches. However this "full" I-S equation is still not the most general 2nd order theory, and in 
a variety of approaches, discussed in Chap. 2, additional 2nd order terms (a vorticity term, a visco-elastic 
term, and more |2031I204[|221| ). each with its own coefficient that takes different values in different theories 
(see Table 2.1), arise. To find the optimal I-S equations for numerical implementation, one needs systematic 
comparisons for these different I-S equations from different approaches. In this section, we will give a detailed 
numerically comparison between the "simplified" and "full" I-S equations [127] . and then briefly study the 
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effects from other higher order terms. 



7.3.1 "Simplified" and "full" I-S equations 
Evolution of momentum anisotropies 

Figure 17.21 shows the temporal evolution of the total momentum anisotropy Sp = ^^^-g^+Tfy) averaged 
over the transverse planc0 for two collision systems (Cu+Cu at 6 = 7fm on the left, Au+Au at 6 = 7fm on 
the right) and three equations of state (EOS I (top), SM-EOS Q (middle), and EOS L (bottom)). The blue 
dashed lines at the top indicate the result from ideal fluid dynamics, the black and orange lines below show 
viscous fluid dynamical results. The black lines show solutions of the "full" I-S equation, the orange ones 
for the simplified TS approach; in each case several values of the kinetic relaxation time t^^ are explored. 
Note that our "full" I-S equation used here is defined by eq. (|3.3p (already written in 2+1 dimensions), but 
not the full I-S equation defined by eq. (|2.3ip in Chap. 3, which strictly holds only for conformal fluids. We 
have, however, tested the two expressions on the left and right side of eq. (|2.32p against each other also for 
the other two equations of state (SM-EOS Q and EOS L) which are not conformally invariant, and found no 
discernible differences. Only for a very long relaxation time = 12ri/sT (not shown in Fig. 17. 2p did we see 
for EOS L a difference larger than the line width, with our result for Sp lying slightly above the one obtained 
with the conformal approximation eq. (|2.32p . 

Comparison of the black and orange lines in Fig. 17.21 shows that the sensitivity of the momentum 
anisotropy Sp to the relaxation time is significantly larger for the simplified I-S equation (orange) than 
for the "full" I-S equation (black), and that the TTr-dependence of Sp even has the opposite sign for the two 
sets of equations. With the "full" I-S equations, Sp moves slowly towards the ideal fluid limit as increases 
whereas with the simplified I-S equations Sp moves away from the ideal fluid limit, at a more rapid rate, 
resulting in a larger viscous suppression of the momentum anisotropy. In the limit r^r — )■ 0, both formulations 
approach the same Navier-Stokes limit. The difference between "full" and simplified I-S equations is largest 
for EOS I which is the stiffest of the three studied equations of state, causing the most rapid expansion of 
the fireball. For this EOS, the simplified I-S equations allow for the largest excursions of tt™" away from its 
Navier-Stokes limit, causing a significant and strongly TTr-dependent increase of all viscous effects, including 
the suppression of the momentum anisotropy (Fig. 17. 2p and elliptic flow and the amount of viscous entropy 
production (see below). 

For the other two equations of state, SM-EOS Q and EOS L, the difference between "full" and "simplified" 
I-S dynamics is much smaller, ranging from ~ 5% for Au+Au to ~ 15% for Cu+Cu for the largest r^r value 
of 67] /sT studied here. Note that the viscous suppression of Sp is much stronger for the smaller Cu+Cu 
collision system than for Au+Au. For SM-EOS Q and EOS L (which yield rather similar results for Sp, 
with differences not exceeding ~ 10%), the results from the "full" I-S equations (black lines) are almost 
completely independent of r^r, even for the small Cu+Cu system. 

The insets in the two upper panels of Fig. 17.21 illustrate the different TTf-dependences for Ep in the "full" 
and "simplified" I-S formulations, by plotting the value of Ep for EOS I at a fixed time t — tq = 4fm/c 
as a function of t^^. One sees that, for the investigated range of relaxation times, the T,r-dependence is 
linear, but that the slope has different signs for the "full" and "simplified" I-S equations and is much smaller 
for the "full" I-S system. Even though VISH2+1 cannot be run for much smaller values, due to nu- 
merical instabilities that develop as the Navier-Stokes limit = is approached, the lines corresponding 
to the "full" and the "simplified" I-S equations are seen to nicely extrapolate to the same Navier-Stokes 
point, as they should. For SM-EOS Q and EOS L, the corresponding lines may no longer be linear, due 

^"^Note that ep as defined in this section includes the effects from both flow velocity and shear pressure anisotropies J125f . In 

Ref. \125l and Chap 4 denoted it by e' in order to distinguish it from the flow-induced momentum anisotropy ° ^ Jj;;; 

(-'o +-'0 ' 

which is based only on the ideal fluid part of the energy momentum tensor and neglects anisotropies in the local fluid rest frame 
caused by the shear pressure tensor tt""" . In the present chapter we drop the prime for convenience. 
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Figure 7.2: Time evolution of the total momentum anisotropy for two collision systems (left: Cu+Cu; 
right: Au+Au), three equations of state (top: EOS I; middle: SM-EOS Q; bottom: EOS L), and three values 
of the kinetic relaxation time as indicated (dotted, dashed and solid curves, respectively). The insets in 
the two top panels show the r^-dependence of the momentum anisotropy £p at fixed time t — ro = 4 fm/c. 
See text for discussion. 
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Figure 7.3: Viscous suppression of elliptic flow, {v^^^^^ ~ V2^^'^°^^) / v}^°^^ , as a function of {l/S)dNch/dy for 
Au+Au with EOS L, r^r = 'd-q/sT and two values of rj/s as indicated. Solid (dashed) lines correspond to 
using the full (simplified) I-S equations, respectively. 

to phase transition effects, but are still characterized by opposite slopes for the "simplified" and "full" I- 
S approaches, with almost vanishing slope in the full I-S case. This agrees with findings reported in |134[135j . 

Elliptic flovif and viscous suppression of V2. systematics 

In this subsection, we will compare the "simplified" and "full" TS equations by studying the viscous 
V2 suppression with different 77/s as inputs. As already shown in Figs. 15. II and 15.2] the slope of V2IS as a 
function of ^ remains positive for viscous hydrodynamics curves even at the highest collision energies 
(or Co values), i.e. V2/ £ continues to increase and evolve in direction of the asymptotic ideal fiuid limit. This 
implies that at higher collision energies the importance of viscous effects decreases. This observation parallels 
the one made in (125j , namely that with increasing collision energy the pr range increases over which viscous 
hydrodynamic predictions for the single-particle momentum spectra can be trusted. The reason is in both 
cases that with increasing collision energy the time until freeze-out grows, and that (at least for constant 
r}/ s as assumed here and in [125] ) during the later stages of the expansion shear viscous effects are small. 

Figure 17.31 shows this more quantitatively. We plot the fractional decrease of the elliptic flow relative to 
its ideal fluid dynamical value, (wj'^'^'*' — t;2"^'^°"'')/''^2^°'*\ as a function of multiplicity density ^ . Larger 
multiplicity densities lead to smaller viscous suppression effects. Larger viscosity results in stronger suppres- 
sion of the elliptic flow. The suppression effects are weaker if the "full" I-S equations are used than in the 
"simplified" approach of Ref. [125] (which also suffers from strong sensitivity to T^r). For minimal viscosity, 
7//s = I/47r, the pT-integrated elliptic flow V2 in Au+Au collisions at RHIC is suppressed by about 20%. 
The suppression is larger at lower energies but will be less at the LHC. 

Entropy production: 

In this subsection, we compare the entropy production from viscous hydrodynamics using the "simpli- 
fied" and "full" I-S equations. Figure [7^ shows that, when the "simplified" I-S equations are used, entropy 
production depends very sensitively on the kinetic relaxation time r^, approaching the much smaller and 
almost completely T7r-independent entropy production rates of the "full" I-S framework in the limit r^r — » 0. 
The large amount of extra entropy production for non-zero in the "simplified" I-S approach must thus 
be considered as unphysical. This is important because this artificial extra entropy production (caused by 
unphysically large excursions of the viscous shear pressure tensor tt™" away from its Navier-Stokes value 
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Figure 7.4: Sensitivity of the entropy production ratio AS /So shown in Fig. 15.41 to the kinetic relaxation 
time Ttt, for the Au+Au colhsion system with eo = 30GeV/fm'^ (corresponding to a collision energy of 
i/i « 200 ^GeV). The three red curves (upper set) are for the simplified Israel-Stewart equations, the three 
black curves (lower set) for the full I-S equations. Solutions with the full I-S equations produce less entropy 
and show very little sensitivity to t^^. 

TT™" = 277(7'"") manifests itself as additional charged hadron multiplicity in the observed final state. Since 
the final multiplicity is used to normalize the initial energy density cq, this causes a significant distortion of 
the initial conditions corresponding to a given set of experimental data, affecting their physical interpretation. 

We conclude that using the "full" I-S equations is mandatory if one wants to minimize artificial effects of 
shear viscosity on entropy production and elliptic flow in the realistic situation of non-zero kinetic relaxation 
times. (We note that, while the value of for the QGP created at RHIC is not very well known, it can 
obviously not be zero). 

7.4 Code verification and comparison between I-S and O-G equations 

Before the code verification process among different groups began, our code VISII2-I-1 had already passed 
several tests of its numerical accuracy: (i) in the limit of vanishing viscosity, it accurately reproduces re- 
sults obtained with the (2-|-l)-d ideal fluid code AZHYDRO [246]; (ii) for homogeneous transverse density 
distributions {i.e. in the absence of transverse density gradients and transverse flow) and vanishing relax- 
ation time it accurately reproduces the known analytic solution of the relativistic Navier-Stokes equation 
for boost-invariant 1-dimensional longitudinal expansion |211j : (iii) for very short kinetic relaxation times 
our Israel-Stewart code accurately reproduces results from a separately coded (2-|-l)-d relativistic Navier- 
Stokes code, under restrictive conditions where the latter produces numerically stable solutions; and (iv) 
for simple analytically parametrized anisotropic velocity profiles the numerical code correctly computes the 
velocity shear tensor that drives the viscous hydrodynamic effects. The details can be found in Ref |126] 
and Appendix lA. 31 

In the rest of this section, we will report on the recent code verification efforts among different groups 
with independently developed viscous hydrodynamic codes. 
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Preliminary results of the TECHQM collaboration (please refer 
to littps://wiki.bnl.gov/TECHQM/index.plip/). I will release the 
material of this part after the related paper has been made available 
on-line. 



7.5 Conclusions 

In this chapter, we first resolved the questions that arose from several recent publications of viscous 
hydrodynamic calculations (with only shear viscosity) which seemed to yield different results I1251I126II134I 
1135] . After accounting for the physical effects arising from different system sizes, different EOS and different 
versions of the I-S equations, the dramatically different V2 suppressions, ranging from 70% to 20%, are largely 
resolved |127| . For a realistic equation of state that implements a quark- hadron transition (here SM-EOS Q 
and EOS L) it turns out that system size effects play a very important role. At RHIC energies and for a 
realistic EOS, the viscous suppression effects for V2/S in Cu+Cu collisions are almost twice as large as for the 
larger Au+Au collision system. Non-negligible differences in the amount of viscous V2 suppression arise also 
from details in the EOS, with a smooth crossover as implemented in EOS L giving 25-30% less suppression 
than a first-order transition as in SM-EOS Q. Compared to system size effects and EOS uncertainties, the 
differences between "simplified" and "full" I-S theory are relatively small, affecting the viscous V2 suppression 
at the 10% level relative to each other. The later code verification among different groups eliminated the 
possibility of significant numerical errors in the independently developed codes. Very good agreement is 
obtained when testing the hydrodynamic evolution and the calculated spectra. 

We also numerically studied the effects of implementing different second order theories for causal rela- 
tivistic viscous hydrodynamics. We first compared the the "simplified Israel-Stewart equations" eq. (|2.30p 
used in our early articles [125|,I126] and the "full Israel-Stewart equations" eq. (|2.31|) used in our later 
work [127j and in Refs. |1341I135] . For the "simplified" approach, we found a strong sensitivity of physical 
observables to the presently poorly known kinetic relaxation time t-j^ for the viscous shear pressure tensor 
TT™", in contrast to a much weaker and basically negligible r^-dependence in the "full" approach. For non- 
zero T^r the "simplified I-S equations" allow for large excursions of tt™" away from its Navier-Stokes limit 
TT™" — 2770-™". These excursions are artificial and disappear in the Navier-Stokes limit t-^ ^ which can, 
however, not be stably simulated numerically. They cause large viscous suppression effects for the elliptic 
flow and large amounts of extra entropy production (i.e. extra final hadron multiplicity). From our study 
we conclude that the "simplified I-S approach" should be avoided, and that a reliable extraction of i]/s from 
experimental data mandates the use of the "full Israel-Stewart equations" |134lll35lll27j . We also studied 
the effects from other higher order terms beyond the full TS equations, such as the vorticity and visco-elastic 
terms obtained in other approaches. The numerical evidence from our systematic study in Chap. 7.3.2 as 
well as the indications from the code verification among different groups suggest that these higher order 
terms can be neglected in numerical calculations. In short, we recommend to use the "full" Israel-Stewart 
equation (j2.311 12.32[) due to its relatively simple form and its weak dependence on the microscopic relaxation 
time. 
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Chapter 8: Extracting Shear Viscosity from Experimental Data — 

Uncertainty Analysis 



8.1 Introduction 

With the efforts from different groups [HSlfT^fTTniTMlfnSlfraTllTinifTiT] . the elliptic flow has now been 
widely accepted as the key observable to constrain the QGP shear viscosity. With the availability of several 
independently developed causal viscous hydrodynamic codes, we are at the threshold for extracting the QGP 
shear viscosity from experimental data. However, several issues must be clarified before we can do so. These 
include: (1) establishing quantitative uncertainty estimates for the hydrodynamically predicted elliptic flow 
related to uncontrolled uncertainties in the initial conditions for the fireball eccentricity and the remain- 
ing uncertainties in the QCD equation of state (EOS) near and above T,.] (2) inclusion of non-equilibrium 
hadronic chemistry into the EOS below Tc (this is known to affect the distribution of the hydrodynamic 
momentum anisotropy over the various hadron species and thus their V2, but has not yet been included in 
viscous hydrdynamic simulations); (3) exploring uncertainties related to the treatment of the final kinetic 
freeze-out process; and (4) including effects from bulk viscosity, especially near Tc. In this chapter, we 
will try to assess the uncertainties introduced by each of these aspects when extracting the QGP viscosity 
from experimental elliptic fiow data [129j . Resolving these uncertainties will require additional theoretical 
work that goes beyond this thesis. The analysis laid out in this chapter can serve as a guide for such research. 

8.2 Dynamical freeze-out and effects from late hadronic viscosity 

In ideal hydrodynamics, one usually imposes "sudden freeze-out" , i.e. a sudden transition from a ther- 
malized fluid to free-streaming particles, on a hypersurface S(x) of constant temperature or energy den- 
sity [26] . The same algorithm has also been used in most of the existing viscous hydrodynamic calculations 
[134pi25|ll26pi35pi27pi43j . Since viscous hydrodynamics is based on an expansion in small deviations from 
local equilibrium, its validity requires the microscopic relaxation time to be much smaller than the inverse 
macroscopic expansion rate, r^cid-u <C 1. This condition (whose long history is discussed in Refs. |26HI262j 
where it is also applied to ideal hydrodynamics) provides a natural criterium for a dynamical freeze-out al- 
gorithm. Dusling and Teaney [137] implemented it into their viscous hydrodynamics. They find that in this 
case the viscous V2 suppression is dominated by non-equilibrium corrections to the local thermal distribution 
function along the freeze-out surface [137] . This is not a collective effect arising from anisotropics of the flow 
velocity profile, but a reflection of non-equilibrium momentum anisotropics in the local fluid rest frame. In 
contrast, for isothermal freeze-out we find |126j that the viscous V2 suppression is dominated by the viscous 
reduction of the collective fiow anisotropy, while local rest frame momentum anisotropics play a much smaller 
role. This comparison shows that a careful treatment of the hadronic decoupling process will be required for 
the quantitative extraction of rj/s from elliptic fiow data. Dusling also found that dynamical freeze-out can 
increase the slope of the multiplicity dependence of the eccentricity-scaled elliptic fiow V2/e [208j . This is 
an improvement over the scaling behavior found in [127j for viscous hydrodynamics with constant rj/s and 
isothermal freeze-out which features a slope that is too small. 

There are other reasons why a proper kinetic treatment of the late hadronic phase is important. By 
matching a realistic hadron rescattering cascade to an ideal fiuid description of the QGP and hadronization 
stages, Hirano et al. |229j showed that the HRG phase is highly viscous and strongly suppresses any buildup 
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of elliptic flow during the hadronic stage. This is consistent with a recent analysis by Demir and Bass [169] 
who found large shear viscosities for their hadronic UrQMD cascade even close to Tc (between 5-10 times 
above the KSS bound). 

To extract the QGP viscosity from elliptic flow data, one needs a good description of the highly viscous 
hadronic stage in order to properly include the additional viscous V2 suppression from the late hadronic 
stage. However, viscous hydrodynamics with a temperature dependent rj/s (using the results from |169| ) 
can not self-consistently describe the highly viscous hadronic stage. The HRG viscosities found in [169j are 
so large that the viscous hydrodynamics description breaks down even during the early time of the fireball 
evolution [263j as one reaches the skin of the fireball where the matter is in the HRG phase. This requires 
the future development of a hybrid approach that connects viscous hydrodynamics (for the QGP phase) with 
a hadron cascade. Moreover, the correct description of the beam energy and centrality dependence of V2 is 
crucially affected by the changing relative weight of QGP and HRG dynamics in building elliptic flow (in 
central collisions or at higher energies the system spends more time in the QGP phase than in peripheral 
collisions or at low energies). A realistic kinetic simulation of the hadronic phase and its freeze-out thus 
appears to be indispensable. 



8.3 EoS 

Lattice EoS for the QGP phase: 

The EOS is a necessary input in hydrodynamic calculations. Early ideal hydrodynamic simulations 
generally implemented the so-called EOS Q, which connects the QGP phase with chemical equilibrium HRG 
phase through a first order phase transition [Mll2341IMll264j (see Sec. 3.3 for details). Lattice QCD with two 
light quark flavors (u and d) and one heavier quark flavor (s) predicts, however a smooth cross-over transition 
at small net baryon chemical potential [265j . Using a quasi-particle model, Refs. [266[|267j constructed a 
lattice inspired EOS (qp-EOS) and applied it to ideal hydrodynamics simulations. From comparison runs 
with different EOS as inputs (qp-EOS vs. EOS-Q etc.), the order of magnitude of phase transition effects 
on hadron spectra and elliptic flow were investigated. One found that the differential elliptic flow for heavier 
particles (such as protons) at low pT is sensitive to the order of the phase transition. However, whether 
experimental data favor an EOS with the 1st order phase transition or a cross-over EOS depends on the 
speed of sound near the phase transition as well as the chemical compositions in the subsequent HRG phase 
(as can be seen, e.g. from the different conclusions reached in Refs. [266J and |267j . respectively). It is thus 
necessary to construct a more realistic EOS, using the most recent lattice QCD data for the QGP phase 
followed by a HRG phase in partially chemical equilibrium, and use it as standard input for hydrodynamic 
calculations in order to eliminate uncertainties associated with inaccurate modeling of the EOS. Recent 
developments in lattice QCD simulations for the EoS in 2-1-1 flavor QCD with a physical strange quark mass 
and almost physical light quark masses |268|l2T| make analytical or tabulated constructions of a Lattice QCD 
inspired, realistic EOS possible, and several groups are working along this direction ^ 269) El- 

The other motivation for constructing such a realistic EOS comes from the extraction of the QGP vis- 
cosity. As pointed in Ref. [129|,ll27j (and discussed in Chap. 7.2), this is an 0(25%) effect on viscous V2 
suppression, depending on whether one uses a 1st order or smooth cross-over phase transition EOS. This 
translates into a 0(25%) effect on the extracted value of the shear viscosity to entropy ratio rj/s from ex- 
perimental data. Considering that elliptic flow at present seems to be the only variable from which rj/s can 
be extracted with some accuracy, it is desirable to eliminate the uncertainties from the EOS using our best 
available theoretical knowledge based on lattice QCD. 

Hadron resonance gas EOS with partial chemical equilibrium: 

The so-called lattice EOS (EOS L) used in this thesis and in Ref. \127^ is still a crude one based on old lattice QCD data 
with unphysical quark masses \232f and without extrapolation to the continuum limit. 
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Early ideal hydrodynamic calculations [93ll234[[96ll264| . assuming simultaneous chemical and kinetic 
freeze-out, successfully described both the shape of most pT-spectra and the integrated and differential 
elliptic flow data from central to semi-central Au+Au collisions at top RHIC energies. However, it failed 
to reproduce the particle yields (total final multiplicities) and it underpredicted the pion spectra at low 
Pt- One example is the experimental pion/anti-proton ratio, which could not be reproduced in these ideal 
hydrodynamics simulation, but had to be fixed by hand by introducing a normalization factor. 

Systematic studies of the measured particle yields for large number of hadron species within a statistic 
model indicates that chemical freeze-out happens right after the phase transition at Tch = 160 — 175 MeV [Ml 
l99] . This is much higher than the kinetic freeze-out temperature, Tth = 90 — 130MeV, which is obtained 
from blast-wave model fits to the slope of the px spectra [270]. This separation between chemical and 
kinetic freeze-out motivates the construction an EOS for the HRG phase that implements partial chemical 
equilibrium among the hadron resonances as discussed in Chap. 3.3 |2351I2361I237. 238J . In this PCE-EOS, 
effective chemical potentials for different hadron species (that are introduced to maintain the measured 
particle ratios as the temperature decreases from Tch to Tth) change the relation between temperature and 
energy density, reducing the kinetic freeze-out temperature for the same decoupling energy density. Holding 
the energy density at the value used in the early ideal hydrodynamic studies 93,234,96], this leads to 
steeper transverse momentum spectra. The experimental spectra can be re-fitted by adding some additional 
flow (either from the pre-equilibrium stage or introduced by shear viscosity), which flattens the spectra 
and compensates for the lower freeze-out temperature. However, the non-equilibrium chemical composition 
of the hadronic fireball also affects the distribution of the total momentum anisotropy among the various 
hadronic species at freeze-out. Together with the additional radial flow effects, this leads to an almost 25% 
larger pion elliptic flow |237| compared with early ideal hydrodynamic results and with experimental data, 
opening room for flnite viscosity both in the QGP and hadronic phases. It is clearly necessary to include 
such a large effect also in the viscous hydrodynamics simulations before a quantitative extraction of rj/s from 
experimental data can be attempted. 

As explained in Chap. 3.3, the hadron cascade model [47j naturally incorporates the partial chemical 
equilibrium through its combination of elastic, quasi-elastic and inelastic cross sections among different 
hadron species. Since these cross section are finite, it also accounts for kinetic non-equilibrium (or viscous) 
behavior in the HRG phase, by solving the coupled Boltzmann equations through Monte-Carlo simulations. 

Hybrid approaches, constructed by coupling (3+l)-d ideal hydrodynamics to such a hadron cascade model 
just below Tc, show simultaneously good descriptions of the particle yields, the pT-spectra and the differential 
elliptic flow V2{pt) for different hadron species from Au-|-Au collisions at top RHIC energies [38l[37]. Pure 
ideal hydrodynamic simulations with PCE-EOS give good descriptions of particle yields and spectra, but 
over-predict the elliptic flow. Viscous suppression of additional elliptic flow buildup in the hadronic phase 
is necessary to reduce the chemical non-equilibrium enhancement of pion elliptic flow without affecting the 
hadron yields. These two facts strongly indicate that the chemical non-equilibrium V2 enhancement and 
viscous V2 suppression tend to partially balance each other in the hadronic phase at RHIC energies. This 
will no longer be true at LHC energies, where the elliptic flow fully saturates in the QGP phase [271j such 
that viscous V2 suppression in the hadronic stage is no longer significant, while the V2 enhancement caused by 
the non-equilibrium chemistry in the HRG still persists since the latter reshuffles the momentum anisotropy 
generated in the QGP phase among different hadron species in the HRG phase. 

Considering the above factors, viscous hydrodynamics requires implementation of the PCE-EOS for the 
HRG phase as well as matching to a hadron cascade. Using a partially chemically equilibrated HRG EOS 
is necessary to ensure approximate insensitivity of the final results to the exact value of the switching tem- 
perature Tch, where one switches from viscous hydrodynamics to the cascade model [272]. Such a proper 
treatment of the HRG phase is necessary before extracting the QGP viscosity from experimental data. A 
relatively easy procedure is to directly connect viscous hydrodynamics with a hadron cascade right after 
Tc, where the HRG is still in chemical equilibrium. However, the PCE-EOS helps to generate the correct 
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hadron chemistry even at lower temperatures, thus enlarging the temperature window where one can switch 
from hydrodynamics to the hadron cascade description, without changing the hadron yields. This makes 
the implementation of a more realistic dynamical freeze-out (see Chap. 8.2 ) possible. Such research is 
important, but will have to be left for the future. 

8.4 Bulk viscosity 

Currently, most of the existing viscous hydrodynamic calculations concentrate on shear viscosity effects, 
but neglect bulk viscosity for simplicity. Bulk viscous effects were investigated in Chap. 6 (see Ref. |130) ). 
We found that the negative bulk pressure (11 = —(^6 in the N-S limit) reduces the total pressure (which 
effectively softens the EOS near Tc) , decelerates both longitudinal and transverse expansion, and suppresses 
the build up of flow. As a result, the hadron spectra become steeper, and the elliptic flow is suppressed 
in the low-p^ region. Fig. 16.31 in Chap. 6 showed that even the "minimal" bulk viscosity, as constructed 
there, leads to ~ 25% additional viscous V2 suppression if one uses the N-S initialization. This translates 
into a ~ 25% uncertainty for the extracted value of QGP shear viscosity rj/s when comparing experimental 
elliptic flow data with viscous hydrodynamics results. However, there are no well-defined theoretical results 
for the bulk relaxation time tu- Critical slowing down during the phase transition suggests a larger rn near 
Tc, but the exact analytical form and its maximum value are still unknown. Using a Tn with peak value of 
5 fm/c near Tc, we found that bulk viscous V2 suppression is very sensitive to the bulk pressure initialization 
one inputs, which ranges from 2% (or zero initialization) to ^ 10% (for N-S initialization), if minimal bulk 
viscosity as defined in Chap. 3.6 is used. 

Considering the current theoretical uncertainties on bulk viscosity itself (the peak value of C/s near the 
phase transition and the bulk viscosity in the hadronic phase) and our pour theoretical control over the 
initialization for the bulk pressure and the bulk relaxation time (see Chap. 3.5 for details), it is thus unclear 
how much uncertainty bulk viscosity can bring to the extraction of a shear viscosity value of the QGP. 
Clearly, bulk viscous effects must be taken into account and additional theoretical research is needed to 
reduce the uncertainties of its value and relaxation time. Our analysis showed that 77/s and (/s will always 
combine in their effects on viscous V2 suppression. At present, it is not clear how to sperate them, possibly 
by using additional observables. 

8.5 Glauber model vs. CGC initialization 

We now come to what may turn out to be the most serious road block for precision measurements of 
the QGP shear viscosity: our insufficient knowledge of the initial source eccentricity e. It has now been 
known for a while that the Color Glass Condensate (CGC) model, implemented in the initial entropy or 
energy density profile via the fKLN paramctrizations (see |2291l227j and Chap. 3.2 for references), leads to 
30% larger initial source eccentricities than the popular Glauber model. Ideal fluid dynamics transforms 
this larger source eccentricity into ~ 30% larger elliptic flow. Since the extraction of t]/s is based on the 
viscous suppression of V2, obtained by comparing the measured elliptic flow with an ideal (inviscid) fluid 
dynamical benchmark calculation, a 30% uncertainty in this benchmark can translate into a 100% uncertainty 
in the extracted value for 77/s. This was recently shown by Luzum and Romatschke (see Fig. 8 in [135j ). 
This uncertainty trumps most of the other uncertainties discussed above. Worse, since the initial source 
eccentricity depends on details of the shape of the fKLN proflle near the edge of the distribution where the 
gluon saturation momentum scale Qs becomes small and the CGC model reaches its limit of applicability, 
there is little hope that we can eliminate this uncertainty theoretically from flrst principles. 

Based on their analysis of charged hadron elliptic flow data from the STAR experiment, allowing for a 
20% systematic uncertainty of these data, the authors of [135] found an allowed range < rj/s < 0.1 for 
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Glauber and 0.08 < rj/s < 0.2 for CGC initial conditions. This amount to a 0(100%) difference for the 
extracted r]/s between these two initialization models. Since the analysis in |135| did not include a com- 
prehensive investigation of effects caused by permissible variations of the EOS near Tc, by bulk viscosity, or 
by late hadronic viscosity and non-equilibrium chemical composition at freeze-out (see proceeding sections) , 
one should add a significant additional uncertainty band to these ranges. Furthermore, correcting the ex- 
perimental data for event-by-event fluctuations in the initial source eccentricity |273| may bring down the 
measured V2 values even below the range considered in [135j . Still, we agree with Luzum and Romatschke 
that, even when adding all the above effects in magnitude (ignoring the fact that several of them clearly have 
opposite signs), viscous hydrodynamics with ry/s > 5 x (l/47r) would suppress the elliptic flow too much to 
be incompatible with experiment. 



8.6 Numerical viscosity 

To study the effects from shear and bulk viscosity one must ensure that numerical viscosity is under 
control and sufficiently small. Simply speaking, numerical viscosity comes from the discretization of the hy- 
drodynamic equations for numerical calculation. It causes entropy production even in ideal hydrodynamics 
without shocks and can never be fully avoided. To minimize numerical viscosity, the flux-corrected transport 
algorithm SHASTA [214] employed by VISH2-H1 (and by its ideal fluid ancestor AZHYDRO [246]) imple- 
ments an "antidiffusion step" involving a parameter A called "antidiffusion constant" [274j . For a given 
grid spacing, numerical viscosity is maximized by setting A = 0. In standard situations, the default value 
A = i minimizes numerical viscosity effects |274) . With A = i and typical grid spacing Ax = Ay — 0.1 fm, 
At — 0.04 fm/c, AZHYDRO generates only 0.3% additional entropy in central Au-fAu collisions. This is 
negligible when compared with the 0(10%) entropy production by VISH2-I-1 for a fluid with real shear 
viscosity rj/s = l/An. 

By increasing the grid spacing in AZHYDRO and/or changing A, we can explore the effects of numerical 
viscosity on radial and elliptic flow. We flnd that numerical viscosity has little effect on the development of 
radial flow but reduces V2 in very much the same way as does real shear viscosity. Since we gauge the effects 
of rj/s on V2 by comparing results from VISH2+1 for rj/s 7^ to those for 77/s — 0, we should explore how 
much in the latter case V2 is already suppressed by numerical viscosity. We can do this by setting rj/s = 
and reducing the grid spacing until V2 stops changing (i.e. until we have completely removed all numerical 
viscosity effects on W2)- In this way we have ascertained that for our standard grid spacing numerical viscosity 
suppresses the differential elliptic flow V2{pt) by less than 2%. 



8.7 Conclusions 

While the elliptic flow V2 generated in non-central heavy-ion collisions is very sensitive to the shear 
viscosity to entropy ratio rj/s of the QGP, it is also significantly affected by (i) details of the initialization 
of the hydrodynamic evolution, (ii) bulk viscosity and sound speed near the quark-hadron phase transition, 
and (iii) the chemical composition and non-equilibrium kinetics during the late hadronic stage. Not all 
of these effects are presently fully under control. Recent attempts to extract the specific shear viscosity 
rj/s phenomenologically, by comparing experimental elliptic flow data with viscous hydrodynamics, have 
established a robust upper limit [13511127] 



<5x-5-, (8.1) 
QGP 47r 



tantalizingly close to the conjectured KSS bound 51465, but further progress requires elimination of the above 
systematic uncertainties. Since some of these influence the build-up of elliptic flow in opposite directions, it is 
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quite conceivable that the QGP specific viscosity is in fact much closer to the KSS bound (?7/s)|j^gg = I/Att 
than suggested by the upper limit (jS.ip . Ongoing improvements on the theory side should help to reduce 
or eliminate most of the above uncertainties, bringing us closer to a quantitative extraction of rj/s for the 
quark-gluon plasma. The single largest uncertainty, however, is caused by our poor knowledge of the initial 
source eccentricity which varies by about 30% between models. As shown in |135| . this translates into an 
0(100%) uncertainty for r]/s. It seems unlikely that theory can help to eliminate this uncertainty from first 
principles. It thus appears crucial to develop experimental techniques that may help us to pin down the 
initial source eccentricity phenomenologically, with quantitative precision at the percent level. 
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Chapter 9: 77/5 from Other Considerations 

and Extraction Methods 



9.1 Shear viscosity and Knudsen number 

The Knudsen number approach starts with fitting the experimental multiplicity scaling curve for V2 /e as 
a function of i ^ by a simple formula: 



nyaro i 
V2 _ V2 1 



hydro 

—, — ■ (9.1) 

The above empirical formula are phenomenologically obtained from a two-dimensional parton cascade model 
simulations |275p276|.l277j from which one can directly measure the elliptic flow V2 by sampling the particle 
distribution, and also calculate the Knudsen number iij^. Kq is the free parameter in this formula, and 
from two-dimensional parton cascade model simulation runs one finds Kq — 0.7 for a dilute system. 

Extracting the QGP viscosity requires an estimate of the Knudsen number from experimental data first. 
It is known that elliptic flow V2 is developed at an early stage, roughly at a time scale R/cg for an EoS 
with constant speed of sound Cg (277j . At the time t ~ R/cg, the particle density can be estimated by 
crn ^ ■5^' assuming that total particle number is conserved and the transverse size of the system does 
not vary significantly. Then one can connect the Knudsen number with the final multiplicity per unit area 
{l/S)dN/dy as follows [375]: 

1 R (inR adN , ^ 

After fitting the experimental multiplicity scaling curve for V2/S (for Glauber initial eccentricity and 
CGC initial eccentricity, respectively) with eqs.( 19. ip and ( 19. 2p . one finds V2"^'^° /s = 0.30 ± 0.02 and 
a = 4.3 mb for Glauber initialization, and v^'^^° je — 0.22 ± 0.01 and a — 7.6 mb for CGC initialization 
(using = 1/3) [272. 

The shear viscosity to entropy ratio is estimated using the formula for a classical gas of massless particles 
with isotropic differential cross section |278| : 

-=0.316—. (9.3) 

s can 

One finds 77/5 = 0.19 for Glauber initialization and rj/s = 0.11 for CGC initialization, with the inputs 
T = 200 MeV and n — 3.9 

fm-3 ^277pl . Since the particle density evolves with time, eq. (|9.3p shows that 
a constant cross section and constant rj/s throughout the fireball evolution are mutuUy exclusive. The -q/s 
value extracted from this method thus represents at best an average value, taken over the history of the 
medium. Detailed hydrodynamic simulations in our group [231' show that the evolution history (and thus 
this average) depends on the initialization (Glauber vs. CGC) and the EOS. The estimated rj/s values from 
Ref. [277j thus have to be considered with caution. 



'^^ Knudsen number K is defined as the ratio of the mean free path and the system size, K = For a parton cascade model 
simulated within a two-dimensional Rx X Ry box, the system size is estimated by R= + ■^)~'^/'^ ■ The mean free path A 

is inversely proportional to transport cross section a and particle density n, X = with n = (N is the number of 

test particles in a Monte Carlo simulation). 

^^n = 3.9 fm""^ is estimated at the time when V2 is developed ]275j . 
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Figure 9.1: Left panel: elliptic flow f2(|y| < 1) from BAMPS ^55] (Boltzmann Approach of Multi-Parton 
Scattering) using as = 0.3 and 0.6, compared with the PHOBOS ^M. and STAR ^280^ data. Right panel: 
The shear viscosity to entropy density ratio rj/s at the central region during the entire expansion [55j . rj/s 
values are extracted from the simulations at impact parameter 6 = 4, 6.3, and 8.6 fm.The upper band shows 
the results with = 0.3 and the lower band the results with = 0.6 



9.2 Shear viscosity from the parton cascade model 

In contrast to the macroscopic hydrodynamic approach, the QGP fireball evolution can also be micro- 
scopically described by a Parton Cascade Model (PCM), solving the Boltzmann equation for colliding gluons. 
The earlier parton cascade model [52l[5T] that only included 2^2 elastic gluon collisions was not be able 
to achieve thermal equilibrium and to reproduce the large experimentally measured elliptic flow, unless 
unrealistically large scattering cross sections, far beyond the expectations of PQCD on which the parton 
cascade is based, were assumed |281) . Recently, a radiative parton cascade model including radiative 2^3 
processes and their inverse has been developed to study, on a semiclassical level, the dynamics of gluon 
matter produced at RHIC energies f54ll282|,r283ll284ll285j . This approach has shown that some well-known 
strongly coupled QGP phenomena (such as fast thermalization of the initial non-equilibrium system, large 
elliptic flow and a small shear viscosity to entropy ratio) can be achieved within a framework that is based 
on weakly coupled QCD [286]. 

The left panel of Fig. 19.11 shows the elliptic flow V2 , calculated from the parton cascade model with only 
two body interactions (gg — > gg) and with multi-parton interactions (gg gg, gg — > ggg and ggg — *■ gg, 
denoted as BAMpS in the figure), using different values of the coupling constant as input [55j . One 
finds that large elliptic flow V2 , comparable to the experimental data, can be achieved within the perturbative 
QCD description if the radiative collision processes gg <-*■ ggg are included. 

The shear viscosity ij can be extracted from the Parton Cascade Simulation through the following formula 
(using the N-S approximation) [287j : 

^"5" Ei?tr + |,a,(lnA)- ^'-^^ 

Boltzmann Approach of Multi-Parton Scattering 
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Here n is gluon density, E is gluon energy, Vz = ^ is the gluon velocity along the beam direction and A is 
the gluon fugacity, A = n/rieq. ^ i?*' denotes the total transport collision rate, the expression ol which can 
be found in Ref [287| . 

The right panel of Fig. 19.11 shows the shear viscosity to entropy ratio r]/s extracted from BAMPS, for 
different impact parameters and values. One finds that rj/s depends weakly on the gluon density or 
temperature (indicated by the time evolution on the horizontal axis), but depends strongly on a^. For 
as — 0.6, which best fits the experimental elliptic fiow data, rj/s ~ 0.08, a value that is very close to the 
minimal KSS bound l/47r, is reached. The authors stated, however, that rj/s may increase once hadronization 
and a subsequent hadron cascade are included [286j . 

The authors of Ref. |288j derived an improved expression for shear viscosity within the second order 
Israel-Stewart formalism, for a one-dimensionally expanding particle system. One finds that although the 
improved 1-S expression increases the extracted rj/s, when compared to eq. (|9.4p . the difference is within 
20% for as =0.3-0.6 [287] . 



9.3 Shear viscosity and radiative energy loss 

The estimation of the shear viscosity from radiative energy loss is inspired by the kinetic theory expression 
of shear viscosity rj [278], 

ri^Cp{p)X„,fp^C{p)—, (9.5) 

where p is the particle density, (p) is the average particle momentum and Xmfp is the mean free path. Again, 
the last expression replaces the mean free path by Xmfp = ip<^tr)~^ with particle density p and transport 
cross section crtr- In perturbative QCD and QED, the transport cross section crtr is dominated by small 
angle scattering and can be expressed by 

a.. = ^/&ifci^, (9.6) 

where (i?cm) is the center of mass energy, with (i?cm) ~ 18T^ for a thermal medium with temperature T. 

On the other hand, the transport parameter q (which describes the per unit length energy loss of a fast 
parton traveling through the QGP by gluon bremsstrahlung) , has a theoretical definition associated with 
the transverse momentum broadening of a fast parton through interaction with the medium: 



^'fcifci^- (9-7) 



Here is the differential cross section for elastic scattering in the medium, and p is the medium particle 
density. 

Comparing eqs. (j9.5|) and ()9.7|) . one finds that the shear viscosity to the entropy density ratio ij/s can 
be expressed as |289] : 

— w 1.25 — , (for weak coupling). (9-8) 

s q 

Here one set C — 1/3 |210ll211j and used (p) ~ 3T, and s ~ 3.6/3 for a thermal ensemble consisting of 
massless bosons. 

In Ref. [289j the relationship between shear viscosity rj and jet quenching parameter q was formally 
derived within the weakly coupled QGP picture, for two different circumstances: 1) the distribution of 
quasi-particles satisfies a linearized Boltzmann equation with a soft scattering collision term in leading-log 
approximation for a pure gluon gas |156j : 2) the evolution of the quasi-particle distribution function satisfies 
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a Fokker-Planck equation, which does not contain a collision term but includes the effects of a random color 
field. The latter contributes a small anomalous shear viscosity [290U291) . One finds that the coefficient f .25 
found in eq. (|9.8p holds in both of these cases. 

In the above weak coupling scenario (which describes the quark gluon plasma by nearly massless quasi- 
particles), one estimates that 77/s = 0.12 — 0.24 for a lower value of q^'^ = 1 — 2 GeV^/fm, and 77/s = 
0.008 - 0.024 for a higher value of g^'^^ = 10-30 GeVVfm, with T w MQMeV. The uncertainty of 
the estimated 7]/s mainly comes from the different values of q obtained from different schemes, and the 
error bands in each estimate comes from whether one takes q to be proportional to e'*/"^ or to T^. The 
authors of [289j pointed out that for the larger q value, the 77/s ratio derived from( [9?8l) strongly violates the 
KSS Bound, which they interpreted as a sign of break-down of the weakly-interacting quasi-particle picture 
underlying the relation( l9.8p . 

However, eq. (|9.8p is no longer true in the case of strong coupling. After checking the shear viscosity and 
q from strongly coupled N=4 superymmetric Yang-Mills (SYM) theory and for pion gases (where the strong 
coupling limit of QCD is exhibited by confinement), Majumder, MuUer and Wang [289) also found that 

— ^ 1.25 — , (for strong couphng). (9.9) 

s q 

So in the strong coupling limit, -q/ s can not be usefully derived from a measurement of q. Considering 
eas. (|9.8l) and (|9.9|) together, they conclude instead that "The determination of both shear viscosity rj and jet 
quenching parameter q from experimental data would give a quantitative assessment of the strongly coupled 
nature of the quark gluon plasma produced in heavy ion collisions" |289| . 



9.4 Shear viscosity and heavy quark diffusion 

Heavy quarks are unique probes to study the transport properties of the QGP. Based on the property 
that the masses {Ale = 1-3 GeV and Mb = 4.2 GeV) of heavy quarks are much larger than the medium 
temperature (T ^ 300MeV), Moore and Teaney proposed a simple model to study the motion of heavy 
quarks by the following stochastic equations |292j : 

^ = m - VDP., m)Cj{t')) - ks^At " (9.10) 

Here rjo is the drag coefficient, which measures the momentum loss per unit time, and characterizes the 
uncorrelated random momentum kicks from the medium. 

The diffusion coefficient D is defined as the probability of "starting a heavy quark at a:: = at i = and 
finding it at the mean squared position Dt at a later time" : 

{x,{t)xj{t)) ^2DtS,j — > 6Dt^{x^{t)). (9.11) 

One finds that D is related to the drag coefficient rye at zero momentum [292) : 

D= ^— . (9.12) 

After perturbatively calculating the heavy quark diffusion coefficient D and comparing it with the weakly- 
coupled QGP shear viscosity 77 at leading order, Moore and Teaney found that the heavy quark diffusion 
coefficient is related to the hydrodynamic diffusion coefficient by a factor of ~ 6 ^292^ , For later conve- 
nience, we rewrite this as follows (using e+ p = Ts): 

'^f=i\TD (wQGP). (9.13) 
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The heavy quark diffusion coefhcient can also be calculated in strongly coupled = 4 Super- Yang-Mills 
theory, using the AdS/CFT correspondence, which gives D ~ l/(27rr) |298j . Recalling the remarkable result 
7//s = l/47r for a large class of strongly coupled theories that have a gravity dual {N = 4 SYM theory is one 
example among them), one finds the following simple relation between the heavy quark diffusion constant 
and the shear viscosity for a strongly coupled system: 

- w —D(2nT) = -TD (N 4 SYM) . (9.14) 
s At: 2 

Fig. l9.2l shows the experimental quenching ratio Raa (see Chap. 1.13) and V2 for "non- photonic electrons" 
from semi-leptonic decays of heavy flavors (i.e. of hadrons carrying charm and botton quarks) for Au4-Au 
collisions at top RHIC energies, together with the theoretical predictions from different groups. The model 
that describes both Raa and V2 well is from van Hees et al. |295) , leading to an estimate of the heavy quark 
diffusion constant D ~ (4 — 6)(27rr). Using the relations (|9.13p and (I9.14p for weakly and strongly coupled 
systems, respectively, one finds: 

V _ ( (1.33 — 2)/47r for weak coupling, , , 

s [ (4 — 6)/47r for strong coupling. 

A more detailed assessment of estimating rj/s from the heavy quark diffusion constant D with different 
models can be found in the recent review article [299] . to which we refer the interested reader. 



9.5 Summary and comments 

In Fig. 19. 3[ we summarize the estimated rj/s from different methods, together with the lattice QCD 
estimates for a pure gluon plasma. One finds that all of these results fall in a relatively narrow band near 
the KSS bound at l/47r. This strongly indicates that the QGP created at RHIC energies is an almost perfect 
liquid. 
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Figure 9.3: A summary for the estimated rj/s from different method: a) from Knudsen number and mul- 
tiphcity V2 scahng (Ref. [277] and Chap. 9.1); b) from parton cascade model with 2^2 and 2^3 
process (Ref. [55] and Chap. 9.2); c) from transport parameter q of radiative energy loss (Ref. |289j and 
Chap. 9.3); d) from heavy quark diffusion constant D (Ref. ^9J and Chap. 9.4); e) from viscous hydrody- 
namics (Ref. |135[|127j and Chap. 8); f) SU(3) gluconic lattice QCD results from Nakamura and Sakai [162] 
and Meyer [163] . respectively. 



When interpreting the error bands in the figure, one must however pay attention to the assumptions 
and limitations of each method. For example, the narrow purple band is from the estimate based on 
parton cascade model, which simulates the dynamics of the pure SU(3) gluons [5S]. To compare with the 
experimental data, it directly converts partons to hadrons, but totally neglects the effects from the phase 
transition and the succeeding hadronic expansion. Including these, the actual error band from the parton 
cascade model based estimate is very likely much larger than the current one. Although in the future the 
parton cascade model can be connected with the hadron cascade model to account for the highly viscous 
hadronic stage, this will not resolve generic difficulties of cascade models in treating the change of vacuum 
structure that occurs during the phase transition and reproducing a realistic EOS (see Chap. 1.12 for further 
discussion). As a result, the parton cascade model based approach also has difficulties to properly account 
for bulk viscous effects, which, as we showed in Chap. 6, should not be neglected when extracting the QGP 
shear viscosity. 

The Knudsen number approach provides a simple way to estimate the shear viscosity to entropy ratio 
rj/s from the experimental multiplicity scaling data for V2/S. However, it only gives an effective viscosity 
averaged over the whole evolution of the fireball, rather than the QGP viscosity, since the Knudsen number 
formula derived in Ref. [277j only deals with time averaged variables. On the other hand, the framework of 
this approach has also been challenged. For example, questions were raised about the basic assumption of 
this approach, which assumes a constant transport cross section a during the whole fireball evolution and 
then extracts a constant rj/s for a particular collision system. However, the Parton Cascade Model shows 
that a constant cross section a corresponds an rj/s that increases monotonically with time, while a constant 
rj/s corresponds to a monotonically increasing a |300j . This raises questions about the self-consistency of 
the Knudsen number approach. 

The radiative energy loss and heavy quark diffusion approaches investigated both weakly and strongly 
coupled media and showed that the corresponding formulas for rj/s are dramatically different. As emphasized 



88 



in Ref. [289j . instead of using the different formulas to extract ry/s, it is more preferable to use them to 
distinguish the strongly or weakly coupled nature of the QGP medium after both 77/s and q have been 
extracted from experimental data. 

Simulations from current viscous hydrodynamic models show rj/s < 5/47r. This error band is almost the 
largest, when compared with other estimates. However it is a pretty safe results, since it takes into account 
all kinds of uncertainties that currently persist (i.e. initial conditions, EOS, bulk viscosity, viscosity of the 
hadronic stage - see Chap. 8 for details). Compared with other methods, viscous hydrodynamics is a tool 
that directly attacks the problem for extracting the QGP viscosity. In the future it can be systematically 
improved through a controlled analysis of all presently persisting uncertainties. This is expected to further 
narrow down the uncertainty in 77/s by a significant amount. 
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Chapter 10: Summary and Concluding Remarks 



Ideal hydrodynamics has been a great success in describing and (in the case of eUiptic flow) even predicting 
relativistic heavy ion colhsion data at RHIC energies. This has led to the well-known announcement that 
"RHIC scientists serve up the perfect liquid" - the quark-gluon plasma (QGP). To answer the question "How 
perfect is the QGP fluid? "quantitatively, one needs a dynamical framework which allows to include dissipative 
effects accurately and consistently. Causal viscous hydrodynamics is such a tool that not only improves the 
description of heavy-ion fireball dynamics in regions of parameter space where ideal hydrodynamics begins to 
fail, but also allows to extract even small values of rj/s through detailed comparison with precise experimental 
data. 

This thesis focused on causal viscous hydrodynamics for relativistic heavy ion collisions. The work in- 
cluded developing a numerical code to solve the causal viscous hydrodynamic equations and searching for 
possible signals for extracting the QGP viscosity. Based on the explicit form of the general 2nd-order Israel- 
Stewart (I-S) equations in 2-|-l-dimensions (derived by my advisor U. Heinz and myself [124] ). assuming 
longitudinal boost-invariance but arbitrary dynamics in the 2-dimensional transverse plane, I subsequently 
developed a (2-|-l)-d viscous hydrodynamic code, called VISH2-I-1 (for "Viscous Israel-Stewart Hydrodynam- 
ics in 2-1-1 dimensions") |125''126'. This code has not only passed extensive tests against the well-established 
(2+l)-d ideal fluid code AZHYDRO and analytically known solutions for problems with reduced dimension- 
ality during the early stage of my work |126j (see Appendix A. 2), but recently also passed detailed code 
verification tests (within the TECHQM Collaboration) against two other viscous hydrodynamics codes that 
were developed independently at the University of Washington and at Stony Brook later on (see Chap. 7.4). 

Using VISH2+1, we numerically studied the effects from shear viscosity |125|ll26j (Chap. 4) and bulk 
viscosity |130] (Chap. 6) on the hydrodynamic evolution of the QGP flreball, the flnal hadron spectra, 
and their elliptic flow coefficient V2- We found that shear viscosity reduces the longitudinal expansion, 
but accelerates transverse expansion, which leads to larger radial flow and flatter spectra when compared 
with an ideal fluid with identical initial and final conditions. Bulk viscosity reduces both longitudinal and 
transverse expansion, resulting in smaller radial flow and steeper spectra. Both shear and bulk viscosity 
suppress elliptic flow V2 at low px- It turns out that V2 is very sensitive to the QGP shear viscosity, and 
that in heavy ion collisions at RHIC energies even the conjectured lower bound, rj/s — l/47r, leads to a 
large suppression of V2. We explored the scaling behavior of V2 with the initial source eccentricity e^, by 
computing V2/Sx as a function of charged hadron multiplicity in both ideal and viscous hydrodynamics, 
for Cu-f-Cu and Au-|-Au collisions at a variety of impact parameters and collision energies. We found that 
smaller systems and collisions at lower energies feature stronger shear viscosity effects, and that non-zero 
shear viscosity breaks the scaling of V2 1 £x with charged multiplicity by small amounts that can perhaps be 
used to further constrain the QGP shear viscosity experimentally [127] (Chap. 5). 

The first attempt to extract the QGP viscosity from RHIC elliptic flow data using viscous hydrodynamics 
comes from the work of Luzum and Romatschke. Their results indicate |135lll29j 

^ <5x-L. (10.1) 

s QGP 47r 

A more precise extraction of the QGP viscosity, however, requires to at least consider the following aspects 
in the future [129] (Chap. 8): 

• Better understanding of the initialization: the largest contribution to the present uncertainty range for 
r\l s comes from uncontrolled uncertainties in the initial conditions for the flreball eccentricity, more 
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specifically, the Glauber initialization vs. the Color Glass Condensate initialization. This issue can 
not be solved by hydrodynamics itself, but requires new experimental techniques to help us pin down 
the initial source eccentricity phenomenologically. In addition, it is important to implement event-by- 
event fluctuations in the initial source eccentricity in viscous hydrodynamic simulations to account for 
additional contributions to the elliptic flow V2- 

• Constructing a better EOS: The non-equilibrium chemistry in the hadronic phase is known to affect 
the distribution of the hydrodynamic momentum anisotropy over the various hadron species and thus 
their elliptic flow, but has not yet been included in viscous hydrodynamic simulations. With the 
availability of new, more accurate lattice QCD data on the EOS, using almost physical quark masses, 
it is now possible and necessary to construct a new EOS that connects this latest lattice EOS above 
Tc with the non-equilibrium hadronic EOS below Tc- 

• Viscous hydrodynamics -|- hadron cascade: The hadronic stage is highly viscous, to the extent that 
the framework of viscous fluid dynamics may break down. It is therefore necessary to match viscous 
hydrodynamics to a microscopic hadronic cascade below Tc. It is also important to investigate how 
well the late non-equilibrium hadronic stage can be reasonably simulated by viscous hydrodynamics 
with a temperature-dependent viscosity/entropy ratio that increases rapidly as the matter approaches 
freeze-out, and explore possibilities of extracting the viscosity of hadronic matter. 

• Including effects from bulk viscosity: My recent research has shown that bulk viscosity also suppresses 
U2, thus adding to the effects from shear viscosity. It is therefore important to include bulk viscosity 
effects when extracting the QGP shear viscosity. However, relaxation times and initial values for the 
bulk pressure are required inputs in viscous hydrodynamic calculations, in addition to the transport 
coefficients and the EOS. Near T^, the bulk viscosity C, can exceed the shear viscosity r] of the strongly 
interacting matter. If the relaxation time rn for the bulk viscous pressure 11 is short, it quickly loses 
memory of its initial valaue, but the relatively large peak value of C/s near Tg can lead to a significant 
viscous suppression of the elliptic flow V2, competing with shear viscous effects. If rn grows rapidly 
near Tc, due to critical slowing down, the bulk viscous suppression effects on V2 depend crucially on 
the initial value of 11: If 11 is zero initially, bulk viscous effects on V2 are almost negligible; if 11 is 
initially large, however, as for the case of the N-S initialization, it remains relatively large throughout 
the evolution, suppressing the buildup of elliptic flow at a level that again competes with shear viscous 
effects. Additional research on initial conditions and relaxation times for the bulk viscous pressure is 
therefore necessary for a quantitative extraction of tj/s from measured data. 

In short, the field of viscous relativistic fluid dynamics is experiencing a fast and healthy development, 
with contributions from and collaborations among several different groups. As a result, we have reached a 
new threshhold: we are in the process of quantitatively extracting the transport properties of the QGP, a 
new type of matter that we knew almost nothing about even a decade ago. What has been covered in this 
thesis is only the beginning of the story. 
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APPENDIX A: 



A.l Coordinates and Transformations 

The three spatial coordinates and time form a fom' dimensional coordinate system = (t, x, y, z), called 
Cartesian coordinates. Throughout this thesis, we use the metric tensor g^^, = g^'^ = diag(l, — 1, — 1, — 1), 
such that four vectors (x^, for example) transform as follows: 

x^' = {t,x,y,z), Xf, = gf.^x" = {t,-x,-y,-z). (A.l) 

One generally sets the z-axis parallel to the beam direction, and correspondingly calls the {x, y) plane the 
transverse plan (with x pointing in the direction of the impact parameter) . Within the forward light-cone 
\z\ < t, rj — T coordinates x™ = (r, x, y, 77) (with r — \/t'^ — z"^ and 77 = i In f^f ) prove more useful in high en- 
ergy particle and nuclear physics. The metric in this coordinate system reads g™" — diag(l, —1, —1, — l/r^), 
5mn = diag(l, -1, -1, -T^). 

Here we list the transformation between Cartesian and r\ — t coordinates: 




A. 2 Details of the viscous hydro code VISH2+1 

A. 2.1 Expressions for vr"^" and ct'"" 
The expressions for tt™" and tr™" in eq. (|3.3p are 
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Here D = u^dr + u^^c^k + and d ■ u = drU^ + S^jU^ + + — . 



A. 2. 2 Velocity finding 

As shown in [124j . since we evolve all three components tt'^'^, tt"^^, and tt"^^ (one of which is redundant 
due to the constraint 7r'^™u,„ = 0), the flow velocity and energy density can be found from the energy- 
momentum tensor components with the same efficient one-dimensional zero-search algorithm employed in 
ideal hydrodynamics |301[I302] . This is important since this step has to be performed after each time step 
at all spatial grid points in order to evaluate the EOS p(e). 

Using the output from the numerical transport algorithm, one defines the two-dimensional vector M = 
{Mx,My)={T'''=-TT'''',T''y-iT^y). This is (up to the substitution p + H ^ p ) the ideal fluid part of 
the transverse momentum density vector; as such it is parallel to the tranverse flow velocity v±_ = {vx,Vy). 
Introducing further Mq = y^^— tt"^"^, one can write the energy density as 

e = Mo-v^-M = Mq- vj_M, (A.4) 



where v± — \Jv'^+Vy is the transverse flow speed and M = ^ M^+AIy . One sees that solving for e requires 
only the magnitude of v± which is obtained by solving the implicit relation [301[I302|I124] 

_ M 

" Mo+p{e=Mo-v^M) + n' ^ '^^ 

by a one-dimensional zero-search. The flow velocity components are then reconstructed using 

Note that this requires direct numerical propagation of all three components (tt'^'^, tt'^^ and tt"^^) since the 
flow velocity is not known until after the velocity flnding step has been completed. Hence the transversality 
constraint ir'^'^Um = cannot be used to determine, say, tt'^'^ from tt'^^ and tt"^^. However, it can be used 
after the fact to test the numerical accuracy of the transport code. 



A. 3 Tests of the viscous hydro code VISH2+1 

A.3.1 Testing the ideal hydro part of VISH2+1 

When one sets tt"™ = initially and takes the limit 77 = 0, VISH2+1 simulates the evolution of an ideal 
fluid, and its results should agree with those of the well-tested and publicly available (2+l)-dimcnsional 
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Figure A.l: Left: Differential elliptic flow V2{pt) for tt^ from 6 = 4fm Cu+Cu collisions and 6 = 7fm Au+Au 
collisions, using EOS Q. Results from VISH2+1 for 77 = and tt™" = (dashed lines) are compared with the 
ideal fluid code AZHYDRO (solid lines). Right: V2{pt) for 7r~ from Cu+Cu collisions at impact parameters 
6 = 4 and 7fm, comparing VISH2+1 evolution with EOS Q (dashed) and SM-EOS Q (soHd) in the ideal fluid 
limit 77 = 0, 7r"" = 0. 



ideal fluid code AZHYDRO ^246^ . Since VISH2+1 was written independently, using only the flux-corrected 
SHASTA transport algorithm from the AZHYDRO package [24611274 ] in its evolution part, this is a useful 
test of the code. The left panel in Fig. lA.ll shows that, for identical initial and final conditions as described 
in Chap. 3, the two codes indeed produce almost identical results. The small difference in the Au+Au 
system at 6 = 7fm is likely due to the slightly better accuracy of AZHYDRO which, in contrast to VISH2+1, 
invokes an additional timesplitting step in its evolution algorithm. 

When comparing our VISH2+1 results with AZHYDRO we initially found somewhat larger discrepan- 
cies which, however, could be traced back to different versions of the EOS used in the codes (EOS Q in 
AZHYDRO, the smoothed version SM-EOS Q in VISH2+1). In the left panel of Fig. E+J this difference 
has been removed, by running also VISH2+1 with EOS Q. In the right panel we compare VISH2+1 results 
for EOS Q and for SM-EOS Q, showing that even the tiny rounding effects resulting from the smoothing 
procedure used in SM-EOS Q (which renders the EOS slightly stiffer in the mixed phase) lead to differences 
in the elliptic flow for peripheral collisions of small nuclei which exceed the numerical error of the code. 



A. 3. 2 Comparison with analytical results for (0+l)-d boost-invariant viscous 
hydrodynamics 

For boost-invariant longitudinal expansion without transverse flow, the relativistic Navier-Stokes equa- 
tions read [21 Ij (for zero bulk viscosity) 

de e + p + r^TT"" ^ ,^ . 

JT + = 0, A.7 

OT T 

rV"'^-^-. (A.8) 
3 r 
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J _ - - numerical results from (2+l)-d viscous hydro with 
homogenous transverse initial energy densities 
- (0+1 )-d analytic solutions 
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Figure A. 2: Comparison between the analytical temperature evolution for (O-l-l)-d boost-invariant Navier- 
Stokes viscous hydrodynamics (solid line) and numerical results from VISH2-(-l with homogeneous transverse 
initial energy density profiles (dashed line). 



For an ideal gas EOS p= ie~T^ this leads to the following analytic solution for the evolution of the 
temperature |211j : 

t(t) ^ /Toy/3 

To V r / 

To test our code against this analytical result we initialize VISH2+1 with homogeneous tranverse density dis- 
tributions (no transverse pressure gradients and flow) and use the Navier-Stokes identification tt'"" — 2ry(T"™ 
in the hydrodynamic part of the evolution algorithm, sidestepping the part of the code that evolves tt™" ki- 
netically. It turns out that in this case the relativistic Navier-Stokes evolution is numerically stable. Fig. lA.2l 
compares the numerically computed temperature evolution from VISH2-I-1 with the analytic formula (jA.Op . 
for rj/ s = 0.08 and Tq = 360 MeV at tq — 0.6 fm/c. They agree perfectly. 



?lSTr,T, 



(A.9) 



A. 3. 3 Reduction of VISH2+1 to relativistic Navier-Stokes theory for small rj 
and 

Having tested the hydrodynamic part of the evolution algorithm in Appendix IA.3.1[ we would like to 
demonstrate also the accuracy of the kinetic evolution algorithm that evolves the viscous pressure tensor 
components. A straightforward approach would be to take VISH2-I-1, set the relaxation time as close to 
zero as possible, and compare the result with a similar calculation as in Appendix I A . 3 . II where we sidestep 
the kinetic evolution algorithm and instead insert into the hydrodynamic evolution code directly the Navier- 
Stokes identity tt™" = 27717'"". Unfortunately, this naive procedure exposes us to the well-known instability 
and acausality problems of the relativistic Navier-Stokes equations. The suggested procedure only works if a 
set of initial conditions and transport coefficients can be found where these instabilities don't kick in before 
the freeze-out surface has been reached. 

We found that sufficiently stable evolution of the relativistic Navier-Stokes algorithm (i.e. of VISH2-I-1 
with the identification tt™" — 277cr™") can be achieved for standard initial density profiles in Cu-|-Cu collisions 
and the simple ideal gas equation of state EOS I by choosing a very small and temperature dependent specific 
shear viscosity ^ =0.01 2oomcV ~ 2GcV • ^'^^ Israel-Stewart evolution we use a relaxation time which is 
correspondingly short: Ttt = = 0.03 fm/c. 
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Figure A. 3: Differential elliptic flow V2{pt) for gluons from b — 7im Cu+Cu collisions, calculated with ideal 
hydrodynamics (blue dashed line), relativistic Navier-Stokes (NS) hydrodynamics (light blue lines), and 
Israel- Stewart (IS) viscous hydrodynamics with ^ = 2GcV ^^'^ '^^ = 0.03fm/c (red lines), using EOS I. The 
lines for NS and IS viscous hydrodynamics are almost indistinguishable. Solid lines show the full results 
from viscous hydrodynamics, dotted lines neglect viscous corrections to the spectra and take only the flow 
anisotropy effect into account. 

Figure [X73l shows the differential elliptic flow V2{pt) for gluons in 6 = 7fm Cu+Cu collisions evolved with 
these parameters. The dashed line gives the ideal fluid result. The solid and dotted lines show the total 
elliptic flow and the anisotropic flow contribution to V2{pt), respectively, similar to the left panel Fig. 14.91 
There are two solid and dotted lines with different colors, corresponding to Israel-Stewart and Navier-Stokes 
evolution; they are indistinguishable, but clearly different from the ideal fluid result. We conclude that, for 
small shear viscosity 77/5 and in the limit — > 0, the second-order Israel-Stewart algorithm reproduces the 
Navier-Stokes limit and that, therefore, VISII2+1 evolves the kinetic equations for tt™" accurately. 

A. 4 Hydrodynamics vs. blast wave model 

As discussed in Chap. [4731 viscous corrections to the final pion spectra from the hydrodynamic model 
have a different sign (at least in the region pt > 1 GeV) than those originally obtained by Teaney [240j . In 
this Appendix we try to explore the origins of this discrepancy. We will see that the sign and magnitude of 
viscous corrections to the (azimuthally averaged) particle spectra are fragile and depend on details of the 
dynamical evolution and hydrodynamic properties on the freeze-out surface. Fortunately, they same caveat 
does not seem to apply to the viscous corrections to elliptic flow where hydrodynamic and blast wave model 
calculations give qualitatively similar answers. 

Following Teaney 's procedure, we calculate tt"™ in the Navier-Stokes limit tt"™ =2770-"™. We do this 
both in the blast wave model and using the results for ct™" from VISII2+1. For the blast wave model 
we assume like Teaney freeze-out at constant r with a box-like density profile e(r) = edoc^(-Ro~'')j where 
Gdcc = 0.085 GeV/fm'^ is the same freeze-out energy density as in the hydrodynamic model for EOS I, and 
-Ro = 6fm. The velocity profile in the blast wave model is taken to be linear, Ur{r) = aQ-^9{Ro—r), with 
oo = 0.5; freeze-out is assumed to occur at Tdcc = 4.1 fm/c. Rq, uq and rdec are somewhat smaller than in 
Ref. [240| since we study Cu+Cu instead of Au+Au collisions. We concentrate here on a discussion of tt'"'" 
for illustration; the expression for cr'"'" is found in Ref. [124] . Eq. (Allc). While tt'"'' from VISH2+1 differs 



96 



along freeze-out surface 



0.5 



- blast wave model 




: u=(0.5r/R„)e(R„-r) 
- R„=6fm 

>< 1,1, 


, 1 , 



0.5 







- ideal hydro, EOS I 


(bl) 




Cu+Cu (iyyx^^ 




- b=Ofm // / 

- /Y / 




/ / T / 







0.005 



> 

o 



blast wave model 


(a2) 




Tj=4.1fm/c 


, 1,1, 


, 1 


2 4 


6 8 r (fm 



r (fm) 
1 

> 

O A 



r (fm) 



-0.5 



- ideal hydro, EOS I 


(b2) 


71 =2tiv u 


T (fm/c) 
— 1 
2 




4 

6 










- , 1 , 1 , 1 , 1 



r(fm) 



0.5 



0.02 
0.01 

t 

> 

S -0.01 

-0.02 
-0.03 



- ideal hydro, EOS I 


(cl) 


Cu-fCu, b=0fm 






1, 







r(fm) 



- ideal hydro, EOS I 




- Cu-i-Cu, b=Ofm 








7r"=2TiV^'> ^ 




1,1,1 







r(fm) 



Figure A. 4: Top row: Velocity profiles from the blast wave model (left) and from the hydrodynamic model 
with EOS I at fixed times (middle) and along the decoupling surface (right). Bottom row: The corresponding 
profiles for the transverse shear viscous pressure tt'"'' in the Navier-Stokes limit, tt'''' = 277V u'^^ . Calculations 
are for central Cu-fCu collisions, and the curves in the middle panels correspond to the times t = 1, 2, 4, 
and 6 fm/c. See text for discussion. 



from 2rja'^^ due to the finite relaxation time Ttt (see Sec. I4.6p , we have checked that the signs of these two 
quantities are the same on the freeze-out surface so that our discussion provides at least a qualitatively 
correct analysis of the viscous spectra corrections in the two models. 

In Fig. IA.4I we compare the freeze-out profiles for the radial flow velocity and 2T]a^^ from the blast wave 
model. In spite of qualitative similarity of the velocity profiles, the freeze-out profiles of are entirely 

different and even have the opposite sign in the region where most of the hydrodynamic particle production 
occurs (left and right columns in Fig. IA.4[1 . The middle column shows that at fixed times r, the hydrody- 
namic profile for 2r]a^^ shows some similarity with the blast wave model in that 2r](7^^ is positive throughout 
most of the interior of the fireball. What matters for the calculation of the spectra via Eq. (|3.5p . however, 
are the values of 2r]a'^^ on the freeze-out surface E where they are negative, mostly due to radial velocity 
derivatives. This explains the opposite sign of the viscous correction to the spectra in the hydrodynamic 
model and shows that, as far as an estimate of these viscous corrections goes, the blast wave model has 
serious limitations. 



A. 5 Cooling rates for ideal and viscous heavy-ion fireballs 

In the early expansion stage of a heavy-ion collision, shear viscosity leads to a reduction of the longitudinal 
and an increase of the transverse pressure. The reduced longitudinal pressure decreases the work done by 
longitudinal expansion, thereby reducing the initial cooling rate of the fireball. This leads to a somewhat 
increased lifetime of the quark-gluon plasma phase. The increased transverse pressure, on the other hand, 
causes larger transverse acceleration and stronger radial flow of the matter than in ideal fluid dynamics. Due 
to the larger transverse flow, the center of the viscous fireball cools more quickly during the late stages than 
an ideal fluid, thereby slightly reducing the total fireball lifetime until freeze-out for central and near-central 
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Figure A. 5: Central freeze-out times for Au+Au and Cu+Cu collisions with differential initial peak energy 
densities eo = e(0, 0; 5=0) as a function of final charged hadron multiplicity density. Points belonging to one 
series correspond to different impact parameters, with central (peripheral) collisions corresponding to large 
(small) multiplicity densities. Left: ideal ffuid dynamics. Right: viscous fluid dynamics with rj/s = 0.08 and 
Ttt = irj/sT. 



collisions (see Fig. IA.5| ). The viscous fireballs created in peripheral collisions don't live long enough for this 
mechanism to manifest itself; they live longer than their ideal counterparts, due to the decreased initial 
cooling rate arising from the smaller longitudinal pressure (left points in Fig. lA.Sp . 

For an equation of state (EOS) with a first-order quark-hadron phase transition, the mixed phase matter 
is free of pressure gradients and hence not accelerated. Velocity gradients in the mixed phase generate, 
however, viscous pressure components whose gradients continue to accelerate the viscous fluid even while it 
passes through the mixed phase where the thermal pressure gradients vanish. As a result, the viscous fluid 
spends less time in the mixed phase than the ideal one. 

In ideal hydrodynamics, a first-order phase transition with a mixed phase (MP) generates large velocity 
gradients near the QGP-MP and MP-HG (HG — hadron gas) interfaces. In viscous fluids such gradients 
generate viscous pressures that act against building up large velocity gradients. In consequence, all prominent 
structures in ideal hydrodynamics that arise from discontinuities in the speed of sound in ideal hydrodynamics 
are washed out by viscous effects. Shear viscosity thus effectively turns a first-order phase transition into a 
smooth cross-over transition. 
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APPENDIX B: Glossary 



RHIC: Relativistic Heavy Ion Collider 
LHC: Large Hadron Collider 
CGC : Color Glass Condensate 
QGP: the quark gluon plasma 
HRG: hadron resonance gas 
MP: Mixed phase 
EOS : equation of state 

PCE-EOS: partially chemical equilibrium EOS (for the HRG phase) 

N-S equation: Navicr-Stokcs equation 
I-S equation: Israel- Stewart equation 
0-G equation: Ottinger-Gremla equation 

VISH2+1: Viscous Israel-Stewart Hydrodynamics in 2+1 dimensions 
AZHYDRO: AZimuthally asymmetric (ideal) Hydrodynamics in 2+1 dimensions 

5^'': metric tensor 

r: longitudinal proper time 

rj: space time rapidity 

7: Lorentz contraction factor 

u^^: four velocity 

A: atomic number 

b: impact parameter 
R: nuclear radius 

eo(r = 0, 6 = 0): initial peak energy density value at r = 

Cdec- decoupling energy density 

tq: initial time 

Tf. freeze-out time 

Cs'. speed of sound 

Tcr' critical temperature 

Tdec- decoupling temperature 

Ex' spatial eccentricity 

Spi momentum anisotropy 

vt' transverse flow velocity 

Nch- produced charged hadrons 

y: momentum rapidity 

Pt- particle transverse momentum 

V2- elliptic flow coeflacient 

S: overlap area 

T^'^: energy momentum tensor 
A'''': conserved charge flow 
f{x,p): distribution function 
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fo{x,p): equilibrium distribution function 

df{x,p): non-equilibrium part in the distribution function f = fo + 5f 

e: energy density 

p: pressure 

T: temperature 

s,S: entropy density, total entropy 
A<S entropy production 

tt''^: shear pressure tensor 
11: bulk pressure 
g'*: heat flow 

bulk viscosity 
77: shear viscosity 
A: heat conductivity 

Ttt: relaxation time for shear pressure tensor 

Tn: relaxation time for bulk pressure 

Ai, A2, A3: other 2nd order transport coefficients in the viscous equation for tt'"' 
^ij.v^ ^(ij^yv) ■ velocity stress tensor 
0'*'^: vorticity tensor 
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